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CHAPTER  1 


Steady-States  in  Two-Species  Particle  Aggregation 

1.1  Preface 

This  chapter  is  concerned  with  stability  analysis  of  configurations  formed  by  a 
large  number  of  two  species  of  pairwise- interacting  particles  in  M2  [60].  The  crux 
of  our  stability  analysis  comes  down  to  the  aggregation  equation 

ut  +  div(u/c)  =  0, 
v  =  -Vif  *  u, 

which  is  used  to  approximate  the  damped  collective  particle  motion  by  the  dy¬ 
namics  of  a  continuum  mass  u(x )  subject  to  the  potential  K. 

To  illustrate  the  mathematical  tools  used,  we  first  introduce  PDE  limits  for  dis¬ 
crete  systems,  discuss  the  aggregation  equation,  give  an  example  of  the  technique 
in  action,  and  then  proceed  to  the  specifics  of  our  project  in  section  1.2. 

1.1.1  PDE  Limits 

A  main  step  of  this  chapter  is  the  use  of  a  nonlocal  PDE  to  describe  aggregate 
particle  interactions  as  the  number  of  particles  goes  to  infinity.  Analysis  of  the 
system  of  particles  is  then  shifted  from  high-dimensional  coupled  ODEs  to  a  PDE. 
For  our  purposes,  this  continuum  limit  is  partially  an  ansatz  because  we  consider 
potentials  for  which  the  particles  remain  contained  in  a  bounded  set  regardless 
of  their  number.  Therefore,  we  assume  from  the  start  that  the  problem  may 
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be  formulated  in  terms  of  measures;  for  the  case  of  finitely  many  particles,  the 
corresponding  measure  is  a  sum  of  dirac  masses.  Once  the  correct  PDE  limit  has 
been  identified,  it  is  often  easiest  to  think  of  the  case  of  finitely  many  particles 
as  a  special  case  of  the  equations  for  a  general  measure,  i.e.  the  continuum  limit. 
For  more  details,  see  sections  1.1.3  and  1.4. 

Continuum  and  hydrodynamic  limits  appear  in  other  contexts  such  as  random 
interacting  particle  systems,  statistical  mechanics,  and  kinetic  theory.  In  most 
of  these  areas,  identification  of  the  correct  limit  involves  more  complex  scaling 
techniques  and  tools  from  probability  theory.  Discussion  of  these  well-established 
topics  is  outside  the  scope  of  this  dissertation,  but  interested  readers  may  consult 
[13,  48,  74]  and  other  references. 

1.1.2  The  Aggregation  Equation 

The  continuum  limit  of  the  single-species  version  of  the  problem  discussed  in  this 
chapter  is  given  by  the  aggregation  equation  [53]  in  M2: 

ut  +  div(-uu)  =  0, 
v  =  -Vif  *  u, 

which  models  a  scalar  quantity  u(x,  t),  typically  referred  to  as  “mass”,  advected  by 
a  velocity  field  v(x,t).  The  velocity  field  arises  as  the  gradient  of  a  potential  K, 
typically  playing  the  role  of  some  kind  of  gravitation,  convolved  with  u.  Thus,  the 
aggregation  equation  is  commonly  used  to  model  the  overdamped  dynamics  of  a 
quantity  with  nonlocal  attraction/repulsion  specified  by  K.  It  appears  in  models 
of  bacterial  colonies,  swarms,  robotic  control,  and  physical  chemistry.  For  more 
applications  and  models,  see  section  1.2. 

Due  to  the  advective  nature  of  the  aggregation  equation,  its  analysis  (see 
[7,  53,  5]  and  many  others)  has  much  in  common  with  that  of  the  fluid  Euler 
equation  [62],  In  the  sense  of  the  Helmholtz  Decomposition  for  vector  fields,  the 
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Euler  and  aggregation  equations  are  orthogonal:  in  the  Euler  equation  the  velocity 
held  is  divergence-free,  while  in  the  aggregation  equation  it  is  the  gradient  of  a 
potential. 

To  illustrate  this  point,  consider  the  aggregation  equation  with  the  Newtonian 
potential: 

ut  +  div(uu)  =  0, 
v  =  VA_1m  =  ViV  *  u, 


i.e.  v  —  Vd  where  A 0  =  c o.  The  Poisson  equation  for  0  is  solved  by  convolution 
with  the  Newtonian  potential  N(x)  =  A-log|x|.  For  more  about  aggregation  in 
this  context,  see  [6]. 

Now  recall  that  the  vorticity  form  of  the  two-dimensional  Euler  equation  is 

[62] 


ojt  +  div(uw)  =  0 
v  =  V^A-1^  =  V±N  *  to, 

i.e.  v  =  where  Aijj  =  to.  Up  to  the  perpendicular  gradient,  the  two  equations 
are  identical.  For  the  vorticity  equation,  the  convolution  defining  the  velocity  is 
known  as  the  Biot-Savart  law  [62]: 


where 


v(x,t)=  /  K2(x  —  y)io(y,  t)  dy 


Ko  =  V±N  =  — 


1 

(  x2 

Xi  \ 

2vr 

X 

2  ’ 

X 

V 

For  the  two-species  problem  discussed  in  this  chapter,  the  continuum  limit  is 
a  generalization  of  the  aggregation  equation  to  two  quantities  with  different  inter¬ 
actions.  Detailed  knowledge  of  the  usual,  single-quantity  aggregation  equation  is 
not  a  prerequisite  for  our  analysis  here,  but  many  of  the  steps  in  section  1.4  will 
be  familiar  to  readers  with  a  background  in  aggregation  or  other  advective  PDE. 
Sections  1.4-1. 6  closely  follow  [87]. 
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1.1.3  An  Example 


To  introduce  the  continuum  limit  technique  in  our  context,  we  consider  the  limit 
of  a  problem  with  a  single  type  of  particle.  Assume  the  particles  occupy  positions 
x-i . . .  Xjv  G  Md,  and  that  every  pair  of  particles  (i,j)  in  the  absence  of  others  will 
move  to  minimize  the  potential  energy 


^(lx*  -xil) 

where  P(x)  is  a  function  with  a  unique  minimum  for  x  >  0.  With  all  particles 
present,  they  jointly  move  to  minimize  the  total  potential  energy 

E(xlr . . .  ,xjv)  E  (i  i) 

To  approach  the  problem,  we  consider  the  continuum  limit  of  (1.1): 

Ec{u)  =  ^  J  J  P(\x  -  y\)u(x)u(y)  dxdy.  (1.2) 

In  the  above  energy,  u{x)  represents  a  continuum  mass  in  space.  If  u(x)  = 
Ef=i^W,  it  is  easily  verihed  that  (1.2)  reduces  to  (1.1)  with  the  assumption 
that  particles  do  not  interact  with  themselves. 

If  we  assume  further  that  u(x,  t )  evolves  the  a  gradient  flow  of  the  energy  (1.2), 
we  arrive  at  the  aggregation  equation 

ut  +  div('uv)  =  0, 
v  =  —  VP  *  u, 


introduced  above.  See  [4]  for  details. 

Alternatively,  we  can  reach  the  continuum  limit  from  the  discrete  energy  E 
directly.  This  is  the  approach  taken  in  sections  1.3  and  1.4  below. 
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1.2  Introduction 


The  collective  behavior  of  systems  of  interacting  particles  gives  rise  to  emer¬ 
gent  phenomena  in  physics,  biology,  chemistry,  and  other  disciplines.  Models 
of  pairwise-interacting  agents  find  applications  in  the  biological  contexts  of  locust 
swarms  [3,  83,  82],  animal  flocks  [25,  55,  63],  and  bacterial  colonies  [85].  These 
mathematical  approaches  to  swarming  have  also  inspired  algorithms  for  coopera¬ 
tive  control  of  robotic  vehicles  [57].  More  questions  for  nonlocal  particle  systems 
arise  in  physical  chemistry:  the  self-assembly  of  nanoparticles  [23,  44]  and  arrange¬ 
ment  of  ions  into  spheres  [58,  59]  are  just  two  examples.  In  the  physical  contexts 
of  granular  gasses  [71]  and  molecular  dynamics  simulations  of  matter  [42],  particle 
systems  also  have  a  central  role. 

All  of  the  above  models,  however  multifaceted,  share  the  same  footing.  Some 
number  of  particles  interact  with  each  other  pairwise  such  that  any  two  particles 
will  repel  each  other  when  they  are  close  and  attract  when  they  are  far;  typically, 
this  attractive  force  disappears  at  very  long  distances.  These  interactions  can 
generate  rich  steady  states  relevant  to  the  models  in  which  they  arise. 

Consider  the  case  when  the  forces  arise  due  to  a  pairwise  interaction  energy 
£(xi,...,xjv)  =  5^-Pdxi-Xjl) 

*7 4? 

where  x,:  denotes  the  position  in  of  the  ith  particle  and  P{r)  is  the  potential 
energy  between  two  particles.  P(r)  is  usually  a  function  with  a  unique  minimum 
such  that  the  force  on  one  particle  due  to  another,  F(r)  :=  — P'(r),  enjoys  the 
repulsive-attractive  properties  mentioned  above.  In  this  framework,  a  steady  state 
pattern  can  be  understood  as  a  minimizer  of  E. 

We  call  the  potential  E  confining  if  its  minimizing  configurations  xi, . . .  ,  xat 
stay  contained  inside  a  compact  set  as  N  — y  oo.  The  question  of  whether  or 
not  a  given  function  P  will  result  in  a  confining  potential  has  been  addressed  in 


5 


terms  of  the  notion  of  H-stability  in  statistical  mechanics;  see  [32],  For  confin¬ 
ing  potentials,  particles  in  ground  states  may  reside  in  space-filling,  co-dimension 
zero  configurations  or  concentrate  on  co-dimension  one  manifolds.  The  question  of 
which  occurs  is  answered  in  [2]  and  the  problem  of  characterizing  ground  states  (or 
steady  states)  has  been  discussed  in  [50,  51,  79,  87],  and  elsewhere.  Applications 
where  both  co-dimension  zero  and  one  solutions  are  of  importance  include  bacte¬ 
rial  colony  growth  under  stress,  point  vortex  theory,  and  the  Thomson  problem 
[1,  12,  22,  85], 

It  is  a  natural  extension  of  the  above  work  to  consider  the  analogous  problems 
for  two  particle  species;  i.e. ,  when  more  than  one  type  of  particle  is  present  in 
the  interactions.  Two-species  models  are  relevant  for  the  phenomena  observed  in 
[59],  where  two  types  of  macroions  in  solution  will  self- recognize  and  assemble 
into  hollow  spherical  structures.  This  self-recognition  of  particle  species  is  a  ro¬ 
bust  phenomenon  observed  in  many  of  the  numerical  experiments  considered  in 
this  chapter.  Two-species  models  also  find  application  in  large  scale  pedestrian 
movement  [73],  and  the  well-posedness  of  said  models  has  been  considered  in  [27]; 
a  general  treatment  of  well-posedness  for  the  two-species  problem  is  given  in  [29] . 
Other  applications  include  opinion  formation  in  groups  consisting  of  ordinary  indi¬ 
viduals  and  strong  leaders  [33]  and  two-species  group  consensus  [36] .  Two-species 
bacterial  aggregation  driven  by  chemotaxis  and  diffusion  is  another  area  of  active 
research,  where  [54]  employs  a  two-species  model  for  localized  vortex  formation  in 
bacterial  colonies.  Global  existence  and  finite  time  blowup  are  considered  in  [24] 
and  [37],  [90]  treats  the  n-species  problem,  and  [45]  and  [80]  discuss  the  stability 
of  uniform  density  and  homogenous  steady  states. 

Our  numerical  experiments  have  revealed  phenomena  which  did  not  appear  in 
the  single  species  problem.  Particle  species  either  mix  or  segregate  based  on  the 
relative  strengths  of  the  inter-species  and  intra-species  forces,  and  occasionally 
settle  in  domains  with  irregular  boundaries  including  cusps.  Asymmetric  steady 
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states  (which  represent  local  minimizers  of  the  potential)  can  be  observed,  and 
nontrivial  structures  form  when  particles  are  H-stable.  These  and  other  features 
indicate  a  substantial  increase  in  complexity  of  the  two-species  problem  over  the 
single  species  problem. 

In  this  chapter,  our  objective  is  to  characterize  steady  states  formed  in  the 
two-species  aggregation  problem  in  the  absence  of  diffusion.  Inspired  by  physi¬ 
cal  [58]  and  numerical  experiments  exhibiting  steady  states  supported  on  or  near 
co-dimension  one  manifolds,  we  wish  to  determine  the  circumstances  under  which 
these  steady  states  form  and  characterize  their  properties  when  possible.  For  ex¬ 
ample,  the  authors  of  [87]  were  able  to  determine  when  the  steady  states  exhibited 
three-fold  or  five-fold  symmetry,  say. 

One  approach  to  this  problem  would  be  to  work  directly  with  the  system  of 
TV]  +  N2  coupled  ODEs,  where  Ni  is  the  number  of  particles  of  species  I  and  N2  is 
analogous.  The  primary  difficulty  is  that  the  number  of  unknowns  increases  as  the 
number  of  particles  increases.  It  is  still  possible  to  pursue  linear  stability  analysis, 
but  the  number  of  linearly  independent  perturbations  to  consider  grows  with  the 
number  of  particles.  Additionally,  linear  stability  analysis  hinges  upon  Ending  a 
meaningful  basis  of  eigen-perturbations  (or  modes)  of  the  system.  This  reduces  to 
an  eigenvector  problem,  but  meaningful  interpretation  of  the  instabilities  becomes 
difficult.  Moreover,  the  results  will  apply  only  for  a  particular  choice  of  N\  and 
N2,  even  though  from  physical  and  numerical  experienmts  we  expect  that  in  many 
cases  the  nature  of  the  instabilities  will  not  change  after  a  certain  number  of 
particles  has  been  reached.  For  example,  a  mode-three  instability  manifesting  as 
a  triangular  arrangement  of  particles  in  the  steady  state  (as  in  [51])  persists  even 
as  more  particles  are  considered. 

An  alternative  method  is  to  consider  the  limit  as  N\ ,  N2  — >  oo.  For  the  single¬ 
species  problem,  this  was  the  approach  used  successfully  in  [51,  87].  The  method 
of  considering  a  continuum  limit  or  hydrodynamic  limit  such  as  this  has  roots  in 


7 


statistical  mechanics,  kinetic  theory,  and  fluids  [13,  48,  74],  When  the  potential 
governing  the  particle  interactions  is  confining,  the  continuum  limit  is  meaning¬ 
ful.  When  the  potential  is  not  confining  for  the  single-species  case,  the  particles 
arrange  into  a  regular  lattice.  In  the  two-species  case,  however,  the  lattice  struc¬ 
ture  formed  may  have  nontrivial  structure  depending  on  the  relative  interaction 
strengths  between  the  two  types  of  particles.  This  phenomenon  is  explored  nu¬ 
merically  in  section  1.8. 

In  this  part,  we  study  the  continuum  limit  of  the  problem  and  the  stability  of 
a  steady  state  consisting  of  two  concentric  rings  of  constant  density.  The  theory 
put  forth  accurately  predicts  instabilities  observed  in  numerical  experiments  and 
the  breakup  of  ring  solutions  into  fully  two-dimensional  patterns.  The  arguments 
presented  here  are  restricted  to  the  two-dimensional  problem,  but  adapting  the 
theory  of  [87]  could  generalize  the  results  to  higher  dimensions. 

1.3  Problem  Description 

Consider  two  species  of  particles,  type  I  and  type  II,  which  occupy  positions 
xi(t), . . . ,  {t),  yi(t),  ■ .  ■ ,  y N2(t)  in  ®12-  Steady  state  patterns  are  minimizers  of 

the  pairwise  interaction  energy 


£(xi,...,xjVl,y1,...,yjv2) 

JVi  N2  Vi  N2 

Y  -xjI)  +  Y  ^(|y*-y;l)  +  EX>(|x*-yil) 

i=i  j= i 


(1,3) 


i,j= 1 

*? 4? 


i,j=  1 

i¥=j 


Vi 


i,j= 1 

*7 


il  2  lx*  —  xil 


N2  , ,  N\  N2 

Yv2(-\yt-yj\ 2)  +EE^(  ^\xi-yj\ 

k  /  i=i  j= i 


i,j= 1 

*7^1 


Above,  Vi(s)  :=  Pi(\/2s),  i  =  1,2,3,  is  simply  a  change  of  variables  to  simplify 


the  calculations. 


We  are  interested  in  gradient  flow  equations  associated  with  (1.3): 


for  j  =  1, . . .  ,N2.  The  right-hand  sides  of  (1.4)  and  (1.5)  have  been  divided  by 
Ni  as  a  simple  rescaling  of  time,  and  in  the  second  line  the  parameter  fi  :=  N2/N\ 
has  been  introduced.  The  factors  1/N\  and  1  /N2  may  also  be  seen  as  normalizing 
each  species  by  its  total  particle  number  or  ‘mass’,  in  which  case  //  represents  the 
relative  mass  of  species  11  to  species  I. 


In  the  above,  g.i(s)  :=  —^r(s)  for  i  =  1,2,3  give  the  ‘forces’  due  to  the 
potentials.  Note  that  gi(s)  is  the  derivative  of  the  rescaled  potential  V  with 
respect  to  its  argument  s  =  1-r2,  where  r  represents  true  particle  distance.  As 
such,  <7j(s)  represents  the  force  only  with  respect  to  the  rescaled  space  variable 
l-r2.  The  true  physical  force —  the  derivative  of  the  potential  with  respect  to  true 
particle  distances  r  and  not  just  with  respect  to  its  argument  |r2 —  has  magnitude 
r<?i(|r2).  The  difference  is  illustrated  in  figure  1.3. 
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One  can  think  of  the  gradient  flow  either  as  an  approximation  to  overdamped 
second  order  physical  dynamics  or  simply  as  a  means  to  identify  minimizers  of 
the  energy.  In  the  next  section,  we  show  that  for  large  numbers  of  particles  the 
gradient  flow  system  (1.4),  (1.5)  may  be  approximated  by  a  nonlocal  PDE  system 
of  advection  equations  similar  to  the  Birkhoff-Rott  equation  for  vortex  sheets  (c.f. 
[62,  78]),  for  which  linear  stability  is  reduced  to  a  sequence  of  eigenvalue  problems. 
Criteria  for  the  stability  of  each  element  in  a  basis  of  perturbations,  and  for  linear 
well-posedness  of  the  concentric  ring  solution,  are  derived.  Numerical  examples 
are  presented,  which  demonstrate  strong  agreement  with  the  theory  put  forth. 


In  this  work  we  consider  the  following  potentials,  which  have  all  been  consid¬ 
ered  in  the  literature  for  the  single  species  problem  [32,  51,  56,  87]:  the  Morse 
potential 


Vi(s)  =  Crie-^s^  - 


power  law  forces 


g.i(s)  =  sPi  -  sqi, 

and  smoothed  step  discontinuity  forces 

tanh|aj(l  —  y/2s)}  +  bt 
9i{s)  =  - — - - - 

with  steady  states  pictured  in  figure  1.1.  Combinations  of  all  three  of  the  above 
types  (and  others)  are  also  plausible;  for  example,  g\  could  arise  from  a  power 
law,  g-2  from  the  Morse  potential,  and  g%  from  the  tanh  force.  See  figure  1.2. 

1.4  The  Continuum  Limit 

For  the  two-species  case,  we  will  say  that  an  energy  E  such  as  (1.3)  is  confining 
if  its  minimizing  configurations  x1; . . . ,  xWl ,  y1; . . .  ,yjv2  stay  contained  inside  a 
compact  set  as  N{  and  N2  — >  oo.  Under  the  assumption  that  the  energy  E 
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Figure  1.1:  Steady  states  of  the  tanh  potential  with  n  —  1,  a\  =  <22  =  10, 
bi  =  b-2  =  0.1,  a3  =  2  and  63  =  —0.7, —0.5, —0.3, —0.1,  0.1,  and  0.3  (from  top 
left  to  bottom  right).  Each  steady  state  consists  of  1000  white  particles  and  1000 
black,  with  (1.4),  (1.5)  evolved  to  a  final  time  t  =  2000. 
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Figure  1.2:  Steady  states  with  g\  and  r/2  tanh  forces  and  a  Morse  force.  On  the 
left,  a\  =  a-2  =  10,  61  =  62  =  —0.3,  Ca  =  la  =  Cr  =  1  and  lr  =  0.1.  On  the  right, 
ai  =  02  =  10,  61  =  62  =  0.1,  Ca  =  la  =  Cr  =  1  and  lr  =  0.02. 


Hyperbolic  tangent  force  Morse  force 


Figure  1.3:  Hyperbolic  tangent  (tanh)  and  Morse  forces  r?(|r2)  with  respect  to 
rescaled  space,  and  rgQr2)  with  respect  to  physical  space.  Tanh  parameters 
a  =  10,  6  =  0.1;  Morse  parameters  Ca  =  la  =  Cr  =  1,  lr  =  0.1. 


12 


is  confining,  the  configurations  of  discrete  particles  approach  continuum  spatial 
densities  p\  and  p2. 


As  we  are  interested  in  the  stability  of  concentric  ring  solutions,  we  seek  den¬ 
sities  which  are  supported  along  one- dimensional  curves  T^f)  =  th^cqf)  and 
UM  =  *2(  a,t)  (parameterized  by  a  G  D  C  M)  which  evolve  with  velocity  fields 
V!  and  v2;  that  is, 

<9<l>i  .  .  . 

=  v1($1(a,  t),t) 

=  v2($2(a,t),t).  (1.6) 


The  velocity  fields  Vi  and  v2  are  determined  from  the  respective  densities  by  the 
continuum  limits  of  equations  (1.4)  and  (1.5):  for  x  G  M2, 


vi(x,t)  = 
v2(x,f)  = 


.9i  1  2|X_y 


.92  1 2  ix  —  y 


2  '  (x  -  y)pi(y,t)  +  93^|x  -  y|2  )  (x  —  y)p2(y,  t)  dy 
2 ^  (x-y)p2(y,t)  +  93^lx-y|2 )  (*-y)pi{y,t)dy, 


(1.7) 


where  we  must  assume  that  N2/Nl  — »  p  as  Ni,N2  — >  oo.  The  parameter  p  is 
absorbed  into  p2  and  still  represents  the  relative  mass 

_  Ir2  P2 

H*  r 

Jr2  P i 

To  determine  the  dynamics  of  pi,  p2,  <E*i,  and  <h2  completely,  we  impose  conser¬ 
vation  of  mass: 


+  V-(piVi)  =  0 

^  +  V-(p2v2)  =  0,  (1.8) 

which  is  implicit  in  the  particle  formulation  of  the  problem.  Equations  (1.6),  (1.7), 
and  (1.8)  specify  a  nonlocal  coupled  advection  system  determining  <&!  and  <f>2, 
from  which  p\  and  p2  will  be  recovered  later;  see,  for  example,  [78]  and  [87]. 
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It  is  worth  pointing  out  here  that  the  p^  are  densities  of  measures  singular  with 
respect  to  Lebesgue  measure  on  M2.  Therefore,  they  should  solve  (1.8)  weakly,  or 
in  the  sense  of  distributions.  We  will  assume  there  exist  /*  locally  integrable  on 
M  such  that  for  Borel  sets  E  C  M2, 


Pi(x,t)  dx  = 


fi(a,t )  da 


(a: 


d&i 


da 


da 


{a:  &i(a,t)eE} 

where  p\  admits  the  natural  interpretation  of  the  density  along  the  surface  T,.  It 
then  follows  that  for  ij;  G  (^“(M2  x  [0,  oo] ) , 


■0(x,  t)pi(x.,  t )  dxdt  = 


'o  JD 

poo  (' 

>0  J  D 


if)(&i(a,  t),t)fi(a,  t )  dadt 

1,9$ 

i/>(&i(ai,t),t)pl(a,t) 


da 


dadt. 


One  can  now  integrate  by  parts  from  (1.8)  to  dehne  what  it  means  for  p^  to  be  a 
solution:  for  all  -0  €  C'~(M2  x  [0,  oo]), 


/  o  J  D 


dib 

dadt  =  0. 


Noting  that  (~^  +  v*‘ V?9)  ($,:(«, t),t)  =  one  can  integrate  by 

parts  to  get 


0  =  J 
=  0 


0  JD 

POO 


^(&i(a,t),t))fi(a,t)  dadt 

d 

i/;(&i(a,t),t))—fi(a,t)  dadt, 


Jo  JD 

where  the  boundary  term  drops  out  because  ip  is  compactly  supported.  It  follows 
that  f(a,t)  =  f(a,  0)  =:  f°(a). 


We  also  rewrite  (1.7)  in  terms  of  integrals  along  T j  and  T2: 


id 


vi(x,t)=  /  gi[  ^|x-$!(a,t)|2  )  (x-$i(a,t))/1°(a)  + 


.9.3  (  ^|x-  $2(<M)|2  )  (x-  $2(a,t))/2°(a)  da, 
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v2(x,t)  =  ^c/2^|x-$2(a,t)|2^  (x-  $2(a,t))/2°(a)  + 

^3(^|x-  $i(a,t)|L>)  (x-  $i(a,t))/i°(a)da. 

Appealing  to  (1.6)  then  yields 

d<&i ,  .  .  _  .  ,  . 

=  vi($i  (a,t),t) 

=  I  21Q|*i(M)-M^)|2)  ($i(M)  -  '$1(0',  *))./? (a') 

+  5,3Q|^i(o;^)  -  <£2(a/,t)|2^  ($i(a,i)  -  $2(a', t))/2(a')  da' 

(1.9) 

and 

<9$o 

—  (a,t)  =  v2($2(a,t),t) 

=  J  02Q|$2(<M)  -  $2(0',  t)|2^  ($2(a,t)  -  $2(a/,t))/2°(a') 

+  03Q|$2(M)  -  $i(a',t)|2j  ($2(a,i)  -  $i(a',  t))/i0(a')  da'. 

(1.10) 


The  two  equations  above  determine  With  these  in  hand,  all  that  is  left  is  to 
determine  py  for  this, 


0  =  =  pS^a,t^ 


d$t 


da 


implies  that 


which  is  enough. 


dp! 

dt 


d_  d®± 

s  dt  da 

Pi  d^  > 

da 


Note  that  when  px(x,  t)  =  YaIi  ^(x  ~  xiW)  and  p2(x,  t)  =  ~  Yi(t)), 

the  equations  (1.7)  evaluated  for  V!  at  Xj  and  v2  at  yt  reproduce  (1.4)  and  (1.5) 
exactly,  up  to  the  time  scaling  introduced. 


Before  we  proceed  to  the  next  section  and  linearize  around  the  concentric  ring 
steady  state,  it  is  worthwhile  to  consider  the  existence  of  such  a  steady  state. 
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If  $1  and  $2  parameterize  concentric  circles  of  radii  R\  and  R2  ■  we  can  take 
D  =  [ — 7r,  7r)  and  <&j(s)  =  ©(s)-Rje  1  (as  in  [87])  with 


©(*') 


cos  s'  —  sin  s' 
sin  s'  cos  s' 


With  no  motion  in  time,  equations  (1.9)  and  (1.10)  give 

0  =  [’  9i  Ql*iW  -  4>i(s')|2)  (*i  W  -  <Ms'))/iV)+ 
SaQlSifs)  -  *2<s')|2)  (*iM  -  *2 (s'))/2V)  *'■ 


0 


f_  92  Q  l®2(»)  -  *2(/)l2)  (*2(0)  -  *2(J))flW)+ 
9s(1|*2(»)  -  *i(o')|2)  (*2(0)  -  #,(S'))/lV)*>'- 


The  (constant)  densities  and  radii  must  satisfy 


p\Ri  =  mi 


where  rrit  =  fR2  is  the  total  mass  of  species  i.  In  the  discrete  case,  m*  =  Nt  and 
so  m2/m\  —  N2/Nx  —  //.  Note  that 


f?W)  =  Pi 


d  s' 


=  PtRi 


so  p  =  /2//1,  and  the  above  equations  can  be  rewritten 


0 


[  9i -  *.(PI2)  (*iW  -  *1(0'))+ 

99sQ|*i(s)  -  *2(s')t)  (*iW  -  *2(0'))*'. 


0 


f  W2QIS2W  -  *2(s')|2)  (*2<»)  -  *2<s'))  + 
9sQ|*2(s)  -  *i(s')|2)  (*aW  -  *i(s'))*'. 
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From  the  definition  of 

$i(s)  -  ^j(s')  =  G(s)[RiI  -  RjQ(s'  -  s)}  e1  =  0(s) 
and  so 

|^i(s)  —  ^(s')!2  =  -R2  +  -R2  —  2  RiRj  cos  (s'  —  s). 


R,  —  Rj  cos  (s'  —  s ) 

— Ri  sin  (s'  —  s) 


We  may  cancel  @(s)  from  both  equations  and  reparameterize  the  integrals  so  that 
s  disappears  as  well,  to  reach 


0  = 


/»7T 

/  0l(-R?(l  -cos  s')) 

R\  —  R\  cos  s' 

!  —  7T 

—Ri  sin  s' 

+ 


fRj  +  Rl  D  D  A 

"  R\R2  cos  s 

Ri  —  R2  cos  s' 

V  2  ) 

—R2  sins' 

ds', 


0  = 


n  7T 

/  g,g2(R22(l- coss')) 

R2  —  R-2  cos  s' 

J —n 

—R2  sin  s' 

+ 


(R\  +  Rl  DD  A 

R\R2coss 

R2  —  Ri  cos  s' 

V  2  J 

— Ri  sins' 

ds'. 


The  second  component  of  each  integral  cancels  because  it  is  odd  on  (— 7r,7r),  and 
so  we  are  left  with 


0  =  /  Rigi(R\(l  —  coss'))(l  —  coss')  + 


(1.11a) 


(R\  +  Rl 


2 


—  RfiR2  cos  s' )  (Ri  —  R2  cos  s')  ds' 


0  =  j  nR2g2(Rl(l  -  cos  s'))  (1  -  coss^T  (1.11b) 

RiR2  coss'^  (R2  —  Ri  cos  s')  ds', 


R\  +  Rl 


which  determine  R\  and  R2.  So  long  as  (1.11a)  and  (1.11b)  have  solutions,  the 
concentric  ring  steady  state  exists.  Of  course,  the  integrands  of  (1.11a)  and  (1.11b) 
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must  be  in  Ll{— 7r,7r],  which  is  true  if  gi(t)t1^2  G  L1[0, 1]  and  the  gt  have  no 
singularities  away  from  the  origin.  Assuming  also  that  E(f)f-1/2  e  Z/1  [0, 1] ,  one 
can  define 

F(RuR2)  ■■=  f  -  cos  s'))  +  yh2(^(l  -  cos  s'))  + 

\xV3  (^±±^  -  RtR2  cos  .s'j  d,s’ 

and  observe  that  (1.11a)  arises  as  =  0  and  (1.11b)  as  =  0.  Then  if  gi,  g2, 
and  g3  are  continuous  (except  perhaps  at  the  origin —  for  commonly  encountered 
potentials  such  as  the  Morse  or  Lennard- Jones  potentials,  this  is  the  case)  (1.11a) 
and  (1.11b)  will  be  satisfied  at  a  maximum  or  minimum  of  F.  To  show,  then,  that 
solutions  Ri  and  R2  exist,  it  suffices  to  show  that  F  attains  a  minimum  for  some 

Ri,  R-2  >  0. 

A  general  proof  of  this  fact  is  difficult  because  the  potentials  V)  will  vary,  and 
in  some  cases  a  concentric  ring  solution  may  not  exist.  However,  for  all  cases 
pursued  below,  (1.11a)  and  (1.11b)  have  solutions  R±,  R2  which  do  give  rise  to  a 
steady  state  solution  to  (1.9)  and  (1.10). 

1.5  Linearization  &  Eigenvalue  Problem 

Recall  that  the  rings  have  been  parameterized  as  <E(<s)  =  @(s)-R,;e i  (where  0  is  a 
rotation  matrix).  Consider  now  a  small  perturbation  of  each  ring  in  the  form 

^^E(s)  =  0i(s)(i?jei  +  eAtej(s))  =  <J?j(s)  +  0(s)eAt6j(s) 

so  that  (defining  A*,  B.;,  C,  D.  E.  F) 

6&i(s)  -  6&i(s')  =  ($i(s)  -  $*(s'))  +  eAt[0(s)e,(S)  -  0(s')e,(s')]  =:  A,  +  e^B*, 

5$i(s)  -  6&2(s')  =  ($i (s)  -  $2(s'))  +  ext[Q(s)e1(s)  -  0(s')e2(s')]  =:  C  +  eAtD, 
J$2(s)  -  J$i(s')  =  ($2(s)  -  ^(s'))  +  eAf[0(s)e2(s)  -  0(s')ei (s')]  =:  E  +  eA*F. 
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Linearizing  (1.9)  and  (1.10),  then  canceling  the  factor  of  ext  appearing  in  each 
term  gives 

A0(s)ei(s)  —  J  gi  ^-|  Ax|2^  Bi  +  -Jj-  ^-|Ax|2^  (A1-B1)A1+ 
»(A|C|2)  D  +  ^Q|C|2j  (C-D)Ccta', 

A©(s)e2(s)  —  J  ficj2  (^|A2|~)  B2  +  |A-2|2i  (A2‘B2)A2+ 

92GiE0F+f  G|E|2)(EF,E<is' 

To  simplify  calculations  below,  define  M  =  M{s,s')  :=  0~1(s)0(s/),  and 

u i  =  ©  1  (s)  A,  Uj  =  (/  -  MT)Rie1 

v  =  0~1(s)C  v  =  {RJ  -  i?iA/T)ei 

w  =  0“1(s)E  w  =  {RJ  -  R2M)e1. 

Then 

@-1(s)Aj  —  {I  —  M)  Rid  =  u i  0_1(s)Bj  =  Ci(s)  -  Metis') 

0"1(s)C  =  {RJ  -  R2M)eL  =  v  0“1(s)D  =  Cl(s)  -  Me2{s') 

0"1(s)E  =  {R2I  -  R1M)e1  =  w  0“1(s)F  =  e2(s)  -  Me^s') 

and  because  0  and  M  are  unitary, 

Aj-Bj  =  [0(s)(7  -  M)Rie1]  .[©(s)ei(s)]  -  [0(s)(J  -  M)Rle1]-[Q{s)Mei{s')} 

=  (/  -  M)Rlerel(s)  -{I-  M)RierMei{s') 

=  (/  —  A/I)Rie\  ■ ei{s )  +  (/  —  MT)Rie1-ei{s'), 

=  Ui-ei(s)  +  ui-ei(s/), 

C  D  =  {RJ  -  R2M)e1-e1{s)  +  {R2I  -  i?iA/T)ere2(s') 

=  v-d(s)  +  v-e2(s/), 

E  F  =  {RJ  -  R1M)ere2{s)  +  {RJ  -  f?2MT)ere1(s') 

=  w-e2(s)  +  w-ei(s'). 
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Finally, 


|Aj|  =  |0  1Ai\  =  |  u.j 

IQi  =  le-'Cii  =  |  Vi  | 

|Ej|  =  |0_1Ej|  =  |wj 


In  terms  of  these  quantities,  multiplying  the  linearized  equations  by  0  1  and 
collecting  terms  multiplied  by  Cj(s)  and  ej(s')  leaves 


Aei(s)  = 


0i 


ui 


M  3 


G1"’) 


r  dg  ifl 

/+M  2|U1 

dg:i  f  1 


/  +  ^,_M2lv0v 


Ui  (8)  Ui+ 

€i(s)  ds' 


(1.12) 


+ 


dg-.i  (  i 


U!  <g)  Ui 


m[  dvl  )  M  + 


El  (s')  ds' 
e2 (s')  ds ', 


Ae2(s)  = 


dgz  f 


lig2  x|u2|  I  +  n—  -|u2|  U2(g)U2  + 


ds  \2 


(1.13) 


dg 3  /: 


+ 

+ 


03  (  „|w|  )  /  +  J  W®W 

11  '2)m  +  ^(1u2|2 


e2(s)  ds' 


-/192  (  2 1  u2 


y 


ds  \2 
d.9.3  / 1 


u2  <g)  u2 


e2(s')  ds' 


03 (  dwl  )  M  +  -^-1  1 w 0 w 


Ci(s')  ds'. 


Explicitly, 


—  |uj(s  -  s')|2  =  H2(  1  -  cos(s  -  s')) 


„/\|2 


-vs— s  =-ws— s 


„'\|2 


—  —  [if2  +  H2  —  2RiR2  COs( 


s  —  s 


Ui  =  Ri 

1  —  cos(s'  —  s) 

,  Ui  =  Ri 

1  —  cos(s'  —  s) 

—  sin(s'  —  s) 

sin(s'  —  s) 

Hi  —  H2  cos(s'  —  s) 

Hi  —  R-2  cos  (s'  —  s) 

V  = 

— H2  sin(s'  —  s) 

,  v  = 

H2  sin(s'  —  s) 
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i?2  —  Rl  COs(s'  —  s) 

R-2  —  Ri  cos(s'  —  s) 

w  = 

—R\  sin(s'  —  s) 

,  w  = 

Ri  sin(s'  —  s) 

and  u  ®  v  denotes  the  matrix  with  i ,  j  entry  UiVj. 

Note  that  all  the  above  matrices  have  even,  periodic  entries  along  the  diagonals 
and  odd,  periodic  entries  off.  With  this  in  mind,  consider  (1.12)  rewritten  as 

p2iz  /*27r  p2ir 

Xei(s)  =  /  M1(s'—s)ds'e1(s)+  M2(s'—s)e1(s')ds'+  M3(s'~ s)e2(s')  ds' 

Jo  Jo  Jo 

(the  superscripts  are  used  to  distinguish  matrices,  not  as  powers)  where  the  di¬ 
agonal  entries  of  M1,  M2,  and  M3  are  even  and  periodic,  and  off-diagonal  entries 
are  odd  and  periodic.  It  follows  that 

p2n 

/  M1  (s' —  s)  ds' ei(s) 

Jo 

is  a  constant  diagonal  matrix  times  €i(s).  For  the  other  two  terms,  we  hope  for 
similar  results  to  yield  an  eigenvalue  problem  in  et  and  A.  Using  the  ansatz 


Xi  cos  (ns) 

yi  cos  (ns) 

el(s)  = 

x2  sin  (ns) 

,  e2(s)  = 

y2  sin(ns) 

similar  to  that  of  [51],  we  compute  the  terms  above  involving  M 2  and  M3: 


s')ei(s')  ds' 


M'n(s  —  s')x i  cos(ns')  +  Mf2(s 
M21(s  —  s')x i  cos(ns')  +  M|2(s 
The  hrst  entry  is  a  linear  combination  of 


s')x  2  sin  (ns')  ds' 
s')x  2  sin  (ns')  ds' 


-2tt 


s)  cos  (ns')  ds  =  /  M21(0)  cos  (n6  +  ns)  d6 


r2i r 


-2t r 


=  cos  (ns)  /  M121(0)  cos(n9)  dQ  —  sin(ns)  /  M21(0)  sin(n0)  dQ 


(*2it 


=  cos  (ns)  /  (9)  cos  (nO)  d6  +  0  (because  Aff,  is  even) 


oc  cos  (ns), 


(1.14a) 


21 


and 


r»27r  r2n 

r2  (  J  n\  J\  j  J  _  /  a  /f 2 


M{2  (s'  —  5)  sinks')  ds'  =  /  Mf2(0)  sin(n0  +  ns)  dd 


f*2n 


r*2n 


=  cos  (ns)  /  M{2(9)  sm(nd)  dd  +  sm(ns)  /  Mf2(0)  cos(nd)  dd 


r»27r 


=  cos(ns)  /  Mf2(0)  sin(n<9)  d<9  +  0  (because  M{2  is  odd) 


oc  cos(ns); 


(1.14b) 


the  second  entry, 


r»27r 


ii4(s' 


s  oc  sin  us 


(1.14c) 


and 


after  similar  calculations.  All  together, 


oc  sin(ns) 


(1.14d) 


2t r 

M2(s  —  s/)ei(s/)  ds1  =  a(n)ei(s) 

where  a  is  a  diagonal  matrix.  The  third  term  will  be  similar  and  will  give  a 
diagonal  matrix  multiple  of  e2,  so  that  the  equation  for  €\  becomes 


A 

Xi  cos  (ns) 

=  a  (n) 

Xi cos  (ns) 

+  b(n) 

yi  cos  (ns) 

X2  sin(ns) 

X2  sin(ns) 

y2  sin  (ns) 

(a  and  b  are  matrix-valued  functions)  and  the  equation  for  e2  is 


yi  cos  (ns) 
i)2  sin(ns) 


c(n) 


yi  cos  (ns) 
i)2  sin(ns) 


+  d(n) 


xi  cos  (ns) 
X2  sin  (ns) 
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Comparing  coefficients  of  cos  (ns)  and  sin(rts)  results  in  an  eigenvalue  problem 
for  xi,x2,yi,  and  y2: 


Xi 

l)\ 

V2 


=  (M1  +  M2(n)) 


X\ 

x2 


+  M3(n) 


M  4(n) 


Xl 

x2 


+  (M5  +  M6(n)) 


Vi 

m 

yi 

V2 


or 


E(n) 


=  A(n) 


with  E  = 


M1  +  M2  M3 
M4  M5  +  M6 


(1.15) 


M1, . . . ,  M6  are  computed  as  in  (1.14a — 1.14d)  and  are  shown  below;  here,  u,, 
v,  and  w  are  functions  of  6.  M1  and  M5  are  diagonal  and  do  not  depend  on  n. 

The  first  two  matrices  determine  the  stability  of  the  species  I  particle  ring  with 
respect  to  frequency  n  perturbations,  with  the  species  II  ring  remaining  fixed: 


Mu  =  /  9i{  ^lui!2)  +^Q|ui|2)(i?i-i?iC°s(s'))2  + 


1 


dg:i  ( 1 


im  (  -|V|JJ  +/U“£i(vhlv|2  )  (^1  -  cos(s'))2  ds' 

M22=  [  ^lQlUl|2)+^QlUl|2)^SinV)  + 


Mn(«) 


M  t2(n) 


M^(n) 


M22(n) 


im (  ^ )  +  d 


l 


=  M 


—  IX 
2 


4/1, 

ds  \  2 

2dgifl 
1  dr  1 

2d9i  ( 1 
i 


Fq  sin2(s/)  ds'. 


-gi  (  —  |u.i |_>  j  cos(6l)  +  R\~i —  (  -|ui|2  )  (1  —  cos(0))2  |  cos (nd)  dd 
2  J  ds  V 2 


gi[  ^lui| 2J  sin(6»)  +  i?f ^2  lUll2  )(1-cos(0))sin(0) 


sin(n$)  d6 


12 


r*7r 


-g i|  2 iUli“ ) cos 


(9)-flilL(i|ui12)  siA°\ 


cos  (nQ)  dO 


=  0 
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after  integrating  by  parts  and  using  (1.11a). 

Due  to  symmetry  in  the  problem,  the  off-diagonal  blocks  M3  and  M1  are 
similar.  M3  represents  the  effect  of  a  perturbation  of  the  species  11  particles  on 
the  ring  of  species  1  particles;  M4,  the  effect  of  a  perturbation  of  the  species  1  ring 
on  the  species  11  ring: 


M3i(n)  =  fj, 
M32  (n)  =  /x 


,9s (  >|2  )  cos(0)  +  %f^|v|2  )  (i?i  -  R2cos(d)) 


ds  \  2 


93^|v|2N)sin(6l)  +  “rfdlv|2)jR2(-Rl“jR2  cos(6>))  sin (9) 


V  2 


ds  \2 


cos (nd)  d6 
sin  [nd)dd 


M|,(n) 


^DlvlA 
ds  \  2  J 


R22  sin2  (6) 


cos  (nd)  d6. 


Mt,(n)  = 
Mj2(n)  = 


-93 1  ;|w|2  |  cos(9)  +  -k  (l|w|2  |  (Ri  -  ft  «*(£>)) 


da  \2 


cos  (nd)  dd 


'  —IT 

r4 


dg3  f 


93  -|w|2  sin(9)  +  —  -|w|- )R1(R2- i?iCOs(6>))  sin(6>) 


M21  (ri)  =  M'(2 


P7T 


M22  (n)  = 


ds  \2 

-.93  (  ^|w|2  )  cos(9)  -  ^Q|w|2  )  R{  sin2(n9) 


sin  (nd)dd 


cos (nd)  dd. 


The  final  two  matrices  M5  and  M6  are  analogous  to  M1  and  M2,  except  they 
determine  the  stability  of  the  species  If  ring: 


Mfi  =  I  ng2(^\u2\2^J  +^Qlu2|2)  (^2  -  R2cos(s'))2+ 
9sQlw|2)  (R2  ~  Ricos(s'))2  ds' 

M22  =  J  /ic/2Q|u2|2^  +9^QlU2l 2j  ^2sin2(s')  + 

*G|W|J) +  ^G|W|J)  - 0 
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after  integrating  by  parts  and  using  (1.11b),  and 


m;2(„) 

MU") 

M«2(n) 


t|u2|2)cos(9)  +  iil^ 
|u2|2)sm(6l)+fl^(t| 


U2 


| U2 1 2 ^  (1  -  cos(O))2 
(1  —  cos(0))  sin(0) 


cos  (nO)  dO 
sin  (nO)dO 


cos (6) 


R*ds 


dg-2  f  1 
\2 


u2| 


Sill' 


\6) 


cos (n6)  dd. 


1.6  Linear  Well-posedness 

Here  we  consider  the  limit  of  the  eigenvalue  problem  (1.15)  as  n  — >  oo.  The  goal 
is  linear  well-posedness ;  that  is,  to  determine  when  the  eigenvalues  A (n)  of  (1.15) 
satisfy  A(n)  <  0  as  n  — >  oo.  That  A (n)  — >  0  as  n  — >  oo  follows  immediately  from 
the  Riemann-Lebesgue  lemma;  the  requirement  that  the  eigenvalues  approach 
zero  from  below  is  important  because  it  demonstrates  that  all  but  finitely  many 
modes  are  stable.  Intuitively,  if  modes  of  arbitrarily  high  frequency  are  unstable, 
the  co-dimension  one  curve  will  break  apart  and  the  density  will  form  a  fully 
two-dimensional  pattern. 


Theorem  1.6.1  (Linear  well-posedness).  Assume  that  the  forces  have  power  se¬ 
ries  representations 


gi(s) 

=  a0spo  +  aispi  +  . . . 

(1.16a) 

02  (s) 

=  b0sgo  +  bisqi  +  . . . 

(1.16b) 

03  (s) 

=  c0sr°  +  cisri  +  . . . 

(1.16c) 

with  ao,bo  >  0  and  cq  ^  0,  valid  in  some  neighborhood  of  the  origin,  where 
Po  <  pi  <  •  •  •  etc.  Define  a  =  and  (3  = 

If  R.\  ^  R-2,  then  the  concentric  ring  solution  to  (1.9,  1.10)  is  linearly  well- 
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posed  if  and  only  if 


a  <  0, 

P  <  0, 


Po  £  2’  U  U  +  1)  2n  +  2), 

9o  £  - ,  0^  U  (2 n  +  1,  2 n  +  2). 


(1.17a) 

(1.17b) 


If  R\  =  R-2,  fj\  =  .92;  and  P  —  1;  7/ie  concentric  ring  solution  to  (1.9,  1.10) 
is  linearly  well-posed  if  and  only  if  (1.17a)  holds  and  either  r0  is  a  nonnegative 
integer  or  r0  >  p0 . 


Before  moving  on  to  the  proof,  a  remark  is  in  order.  Theorem  1.6.1  as  stated 
does  not  cover  the  cases  when  R\  =  i?2  but  g\  ^  <72  or  p,  ^  1.  However,  (1.11a) 
and  (1.11b)  point  out  that  unless  p  =  1  and  g\  =  <72,  it  is  very  unlikely  that 
R\  —  Ra]  for  two  arbitrary  potentials  g\ ,  92  and  ratio  /i,  it  is  a  measure-zero  type 
event  that  the  radii  equations  would  have  such  solutions. 

Proof  of  Theorem  1.6.1.  The  analysis  relies  primarily  on  asymptotic  expressions 
for  the  integrals  occurring  in  M1, . . .  ,  M6,  which  necessitates  the  assumptions  of 
(1.16).  Substituting  (1.16)  into  the  formulas  for  M1, . . . ,  M6  leaves  an  eigenvalue 
problem  where  each  entry  of  E  is  a  potentially  infinite  sum  of  integrals.  However, 
showing  that  E  has  negative  eigenvalues  is  equivalent  to  showing  its  leading  minors 
alternate  sign,  and  it  is  easy  to  see  that  in  each  entry  of  E,  only  those  terms  which 
decay  most  slowly  will  affect  the  eigenvalues  in  the  limit. 

In  practice,  the  values  a  =  Mj^  and  (3  =  must  be  evaluated  analytically 
or  numerically,  because  they  are  independent  of  n.  Other  entries  of  E  may  be 
evaluated  asymptotically,  and  considering  one  of  these  entries  gives  an  idea  of 


26 


how  to  proceed: 


~  a0-RiPo  /  — (1  —  cos  9)Po  cos(9 )  cos (n6)  +  p0(l  —  cos  9)Po+1  cos (n9)  d9 


—  (1  —  cos6))P0-[cos(n  —  1)9  +  cos(n  +  l)d]  + 


=  a0i?iPo 

Po(l  —  cos  9)Po+1  cos (;n9)  d9 


where  we  used  the  trig  identity 


cos(a;)  cos (y)  = 


cos(x  —  y)  +  cos(a;  +  y ) 


By  a  similar  use  similar  identities,  it  turns  out  that  all  entries  of  E  reduce  to  linear 
combinations  of  one  integral: 

/7T 

(c  —  cos  9)p  cos  n9  d,9 , 

-7T 

and  we  are  interested  in  the  behavior  of  /  for  fixed  c  and  p  as  n  — >  oo. 

For  c  =  1  and  p  >  —1/2  (which  is  necessary  for  the  integrals  to  converge),  an 
explicit  formula  with  asymptotics  is  available  from  [87] : 

—C(p)  sin(7rp) 


7(1,  n,p) 


n 


2p+l 


(1.18) 


where  C(p)  >  0  is  a  positive  constant  depending  on  p.  This  asymptotic  form  may 
also  be  arrived  at  by  stationary  phase  analysis. 

For  c  >  1,  it  can  be  shown  readily  via  integration  by  parts  that  I  decays  faster 
than  any  polynomial:  for  any  integer  k,  there  exists  a  constant  C(c,p,k)  such 
that 

\I(c,n,p)\  <  C(c,p,k)n~k.  (1.19) 

For  M2  and  M6,  (1.18)  gives  the  relevant  rates  of  decay.  M3  and  M4  are  more 
complicated  and  the  analysis  breaks  into  cases. 
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Case  I:  R\  ^  R2.  In  this  case  (1.19)  shows  that  M3  and  M4  approach  zero 
faster  than  any  of  the  other  entries  of  E.  and  (1.15)  asymptotically  decouples  into 
two  quasi-single  species  problems 


(M1  +  M2) 

X\ 

=  Kn) 

X\ 

and  (M5  +  M6) 

V\ 

=  A(n) 

X2 

X2 

V2 

V2 

Before  moving  on,  it  is  worth  pointing  out  that  the  decoupling  has  a  pleasant 
physical  interpretation.  M3  represents  the  effects  of  perturbations  of  the  species 
II  particle  ring  on  the  species  I  particle  ring,  and  vice  versa  for  M4.  The  fact  that 
these  vanish  from  (1.15)  as  n  — >■  00  means  that,  when  R\  ^  R2,  high  frequency 
perturbations  of  the  ring  of  species  I  particles  have  no  effect  on  the  eventual 
stability  of  the  ring  of  species  II  (and  vice  versa). 

For  each  of  the  problems  in  (1.20)  to  have  negative  eigenvalues,  it  is  necessary 
and  sufficient  that 


tr(M4  +  M2)  <  0,  det(M4  +  M2)  >  0, 
tr(M5  +  M6)  <  0,  det(M5  +  M6)  >  0. 


(1.21) 


We  turn  first  to  M2:  substituting  in  the  assumed  power  series  representation 
of  g\  and  discarding  all  but  the  most  slowly  decaying  terms, 

/IX 

—  (1  —  cos  9)Po  cos (6)  cos {116)  +  p0(l  —  cos  0)Po+1  cos (nd)  dO 


lll  '  ~  0>qR\ 

=  a0RlPo 

Po(l  —  cos  0)po+1  cos (nO)  dO 
=  a0RlPO 


(1  —  cos6))P0-[cos(n  —  1)6  +  cos(n  +  1)6]  + 


1  1 

-J(l,n  -  l,po)+p0I(l,n,p0  +  1)  -  -J(l,n  +  l,po) 


a0R2^0C{po)  sin(7rp0) 


+ 


L(n  —  1)2po+1  [n  -p  i)2po+i 


a0RPoC{pQ)  sin(7rp0) 

,^2po+l 
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[2 

l12 


71 

0 

0$ 

0 

C3 

PIT 

/  (1“ 

'  —7 r 

a0RlPo 

Po  +  1  7 

2 

Po  +  1 


I(l,n+l,p0) 


a0RlPoC(p0 )  sin(7rpo)(p0  +  1) 


L(n  —  l)2po+1  (n  +  1)2po+i  J 


^  a0R^0C(p0)  sin(7rpo)(p0  +  1)(2 p0  +  1) 

j^2po~\-2 

and  =  M22.  The  final  entry  is  treated  by  integration  by  parts: 


M22  = 


— gi(R\(l  —  cos  9))  cos  9  cos(n9) 


— To  [d i(Ri(l  —  cos0))l  sin#cos(n$)  dO 
du  1 

/IT 

gi(R\{l  —  cos  9))  sin#sin(n$)  d9 


naoR1 

2 

na0RlPO 

2 

?2po 


2po  fir 


/  (1  —  cos#)P0[cos(n  —  1)9  —  cos  (n  +  1)9] 

J  —  7T 

[J(l,  n  -  l,p0)  -  1(1, n  +  1,Po)] 


na0Rf°C(p0)  sin(7r p0) 


L(n-l)2P“+1  (n  +  i)2po+1 


a0RlPOC(p0)  sin(vrpo)(2p0  +  1) 


7j2po+l 


The  analogous  work  for  M6  looks  almost  exactly  the  same,  except  with  qo  replacing 
Po- 

With  those  expansions  in  hand,  one  can  asymptotically  compute  the  terms 
appearing  in  (1.21): 


tr(M1  +  M2)  ~  a 
det(M1  +  M2)  ~  o:M22 


and  so  we  need  only  require  that  a  <  0  and  M22  <  0.  The  asymptotic  expression 
for  the  latter  is  negative  so  long  as  sin(7rp0)  is,  and  this  leads  to  (1.17a).  The 
problem  for  M:)  +  M6  is  nearly  identical,  and  yields  (1.17b)  in  exactly  the  same 
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way.  It  is  worth  noting  here  that  the  criteria  for  linear  well-posedness  are  very 
similar  to  those  for  the  single-species  case  explored  in  [51]  and  [87]. 

Case  II:  R\  =  R-2  R.  As  mentioned  earlier,  this  is  very  unlikely  unless  g\  =  g2 
and  n  —  1;  so,  we  will  assume  that  is  the  case.  Then  M1  =  M5,  M2  =  M6,  and 
M3  =  M4,  so  E  simplifies;  however,  the  rate  of  decay  of  M3  is  not  as  fast  now  and 
so  it  must  be  taken  into  account.  We  determine  when  E  has  negative  eigenval¬ 
ues  by  checking  when  its  leading  minors  alternate  sign.  Asymptotics  for  M3  are 
necessary,  but  M3  has  the  same  form  as  M2  with  g3  replacing  g3  and  R  replacing 
Rv 

M3  c0R2r°C(r 0)  sin(7rr 0) 

11  ~  n2ro+l  ’ 

at3  c0R2roC(r0 )  sin(7rr0)(r0  +  1)(2 r0  +  1) 

M12  =  M21 - - , 

A/t3  c0R2r°C(r 0)  sin(7rr0)(2r0  +  1) 

M22 - - • 

The  first  minor  of  E  is  then  (M1  +  M2(n))n  — >  a  as  n  — >  oo,  so  we  must 
require  that  a  <  0  . 

The  second  minor  is  det(M1  +  M  2M)  rsj  q;M|2  (see  case  I),  so  we  require  that 
(1.17a)  holds. 

The  third  minor  begins  to  include  terms  from  the  cross-particle  interaction 
force  r/3 ,  and  works  out  to  be  (to  leading  order  in  n) 

a2M22  -  a(M32)2  =  — CVn-(2po+1)  +  C2  sin2(7rr0)n“{4ro+4) 

where  C\  and  C2  are  some  positive  constants  with  respect  to  n.  We  have  already 
required  for  the  first  and  second  minors  that  a  <  0  and  M22  <  0,  which  is  why 
the  first  term  is  negative  and  the  second  is  positive.  So  we  must  require  that 

ro  is  a  nonnegative  integer  or  2po  +  1  <  4ro  +  4. 
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The  fourth  minor  works  out  to  be  (again  to  leading  order) 

a2  [(My2  -  (My2]  =  C'in_2(2po+1)  -  C2  sin2  (nr0)n-2{2ro+1) 

where  C\  and  C2  again  denote  positive  constants.  So  we  must  require  that 

r0  is  a  nonnegative  integer  or  r0  >  p0- 

These  restrictions  also  imply  that  the  third  minor  is  negative.  This  last  minor 
gives  the  final  requirement  for  the  double  ring  solution  to  be  linearly  well-posed 
and  yields  the  criterion  in  the  theorem  for  the  Ri  =  R2  case.  □ 

1.7  Numerical  Examples 

All  numerical  solutions  here  and  in  figure  1.1  were  computed  using  a  simple  for¬ 
ward  Euler  scheme  with  an  adaptive  time  step  chosen  as  large  as  possible  while 
requiring  that  the  energy  of  the  system  (1.3)  decreases  at  each  time  step.  A  scheme 
with  higher  order  accuracy  is  not  necessary,  since  we  only  seek  a  minimizer  of  the 
energy  (1.3).  Alternatively,  choosing  a  time  step  based  on  an  estimate  of  the  local 
truncation  error  (as  in  [52])  is  also  efficient  and  yields  the  same  results.  All  initial 
conditions  are  taken  to  be  independently,  uniformly  distributed  on  a  square. 

1.7.1  Mixing  or  Separation  of  the  Two  Species 

The  theoretical  predictions  agree  very  well  with  numerical  observations.  Figure 
1.4  shows  an  example  in  which  the  two  particle  species  may  mix  or  segregate 
based  on  the  relative  strengths  of  the  inter-species  and  intra-species  interactions. 
The  numerical  destabilization  of  the  alternating  ring  structure  and  appearance  of 
mode  two  instability  coincides  exactly  with  the  negative  to  positive  sign  change 
of  an  eigenvalue  of  (1.15)  with  n  =  2. 

Generally,  it  was  observed  that  when  symmetric  intra-species  interactions  g\  = 
g2  are  stronger  near  the  origin  than  the  inter-species  interaction  r/3 ,  mixing  of  the 
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Figure  1.4:  Left:  An  alternating  particle  ring.  Forces  g\ (s)  =  </2(s)  =  1  +  2(1  — s)  + 
.s  — 1  /4  —  0.9357796257,  g3(s)  =  0.5</i(s).  Center:  Separated  particle  ring.  Forces 
r/!  and  r/2  are  the  same  as  on  the  left,  but  #3  =  1.01<7i.  Right:  the  true,  physical 
force  rg\{r2 / 2). 

species  tends  to  occur;  when  g3  is  stronger,  separation  tends  to  occur.  See  hgure 
1.9.  Self-recognition  of  species  (hgure  1.5)  may  be  considered  a  particular  type  of 
separation,  and  inasmuch  as  it  can  be  characterized  by  a  mode  one  instability  it 
can  be  predicted  using  theorem  1.6.1. 

1.7.2  Instabilities  of  Low-frequency  Modes 

Figure  1.5  shows  examples  of  mode  five  and  mode  one  instabilities.  The  alternat¬ 
ing,  symmetric  mode  five  arises  from  interaction  forces  defined  in  terms  of 

G3(s)  =  1  +  (1  —  s)  +  (1  —  s)2,  (1.22) 

G5(s)  =  -(1  -  s)2  +  (1  -  s)3  -  (1  -  s)4, 

G0(s)  =  1  +  2(1  -  s)  +  s~1/A  -  0.9357796257, 
as 

gi(s)  =  G5(s)  +  1.3  •  G0(s), 

<72 (s)  =  G5(s)  +  1.3  •  G0(s), 
g3(s)  =  0.2*G0(s). 
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-1  0  1-2-1012 


Figure  1.5:  Left:  symmetric  mode  five  instability.  The  eigenvector  of  E 
(1.15)  with  positive  eigenvalue  is  [0.02,  —0.71,  —  0.02, 0.71]T;  the  positive-negative 
[a,  b,  —a,  —b]  structure  corresponds  to  the  symmetric  steady  state  observed,  where 
the  species  II  density  is  perturbed  with  the  opposite  sign  as  the  species  I  density. 
Right:  mode  one  instability. 

The  mode  one  instability  is  due  to  the  interaction  forces  defined  as 

gi(s)  =  G0(s), 
g2(s)  =  G0(s), 

g3(s)  =  10 -4(Cre-^lr  -  Cae 

where  GO  is  from  (1.22)  and  g3  Morse  with  Ca  —  1,  la  —  5,  Cr  =  4 ,lr  —  0.5. 

Coupling  effects  of  the  rings  on  each  other  can  be  seen  in  figure  1.6,  which 
shows  coupling  between  type  I  particles  (white)  with  a  mode  three  instability  and 
type  II  particles  (black)  with  mode  five.  The  interaction  forces  are  defined  in 
terms  of  (1.22)  as 

Si (s)  =G3(s)  + 1.1158 -G0(s), 
g2{s)  —  G5(s)  +  1.3  •  G0(s), 
g3(s)  =  - Ks , 
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with  K  =  0,  1,  and  4,  designed  in  [86]  to  exhibit  pure  mode  3  and  5  instabili¬ 
ties.  The  sequential  disappearance  of  the  instabilities  in  figure  1.6  corresponds  to 
the  eigenvalues  of  those  modes  becoming  negative  (i.e.  eigenvalues  from  (1.15)), 
confirmed  numerically. 

Even  when  the  concentric  ring  solution  is  not  linearly  well-posed,  the  stability 
or  instability  of  low  frequency  modes  is  still  borne  out  in  the  ground  state.  Figure 
1.7  shows  the  occurrence  of  a  mode  two  instability  in  such  a  scenario,  which 
is  correctly  predicted  by  the  linear  stability  theory  even  though  the  linearized 
equation  is  ill-posed.  When  multiple  modes  become  unstable,  it  is  possible  that 
one,  several,  or  all  of  the  unstable  modes  will  appear  in  the  ground  state.  The 
question  of  those  which  do,  and  to  what  degree,  is  determined  by  the  particular 
nonlinearity  of  the  problem  and  is  outside  the  scope  of  the  linear  theory.  For 
the  single  species  case,  this  was  pointed  out  in  [86]  and  the  problem  of  which 
modes  appear  is  still  open.  At  this  time,  even  less  is  known  about  the  two-species 
problem.  Weakly  nonlinear  analysis,  considered  in  [79],  may  prove  useful  for  this 
purpose  and  is  one  of  several  considerations  for  future  work. 

1.7.3  Linear  Ill-posedness 

There  are  several  factors  which  limit  the  effectiveness  of  the  linear  theory  when 
the  stable  ground  state  is  far  from  the  concentric  ring  solution.  When  R\  —  i?,2, 
as  is  the  case  in  many  examples  with  g\  —  g2,  the  inter-species  interaction  g3  must 
be  o(s~l/2)  as  s  — »  0  for  the  integrals  in  1.15  to  exist.  This  condition  is  not  met  in 
any  of  the  cases  of  figure  1.1,  and  so  the  theory  may  not  be  applied.  It  is  possible 
to  replace  gi(s)  by  gi((l  +  e)s),  and  if  e  is  sufficiently  small  then  the  observed 
steady  state  is  qualitatively  indistinguishable  from  the  unperturbed  version  but 
Ri  and  i?2  are  no  longer  equal.  The  theory  does  apply  to  this  perturbed  problem, 
but  another  difficulty  arises:  most  or  all  modes  are  unstable. 
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(a)  K  =  0 


(b)  K  =  1 


(c)  K  =  4 


Figure  1.6:  Modes  three  and  five  stabilize  each  other  as  cross-particle  attraction 
increases.  Bottom  right:  true  forces  corresponding  to  g\ ,  £%  and  g%  with  K  =  1. 
Changing  K  scales  the  coupling  force  due  to  #3. 
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Figure  1.7:  Mode  two  becomes  unstable  for  the  power  law  potential.  Left:  g\(s)  = 
g2(s)  =  s~0A5,  g3(s )  =  -s0-15.  Right:  cq(s)  =  g2(s)  =  s~0-15,  g3(s)  =  -s0'2. 

This  type  of  total  instability  also  occurs  in  both  panels  of  figure  1.2.  While  the 
instability  of  all  low  modes  is  certainly  consistent  with  the  observed  steady  states, 
it  is  largely  uninformative.  In  the  majority  of  cases  when  a  mode  1  instability 
appears,  all  low  frequency  modes  are  unstable  as  well  and  the  stable  steady  state  is 
sometimes  asymmetric  (with  the  gradient  flow  coming  to  rest  at  a  local  minimum). 
Nevertheless,  the  linear  stability  theory  may  still  provide  some  insight.  In  the  right 
hand  panel  of  figure  1.2,  each  low  mode  has  one  unstable  eigenvector  except  for 
mode  3,  which  has  two;  however,  only  one  and  not  a  linear  combination  of  both 
appears  in  the  steady  state.  Why  one  and  not  the  other  or  a  combination  of  both 
appears  is  impossible  to  determine  by  the  linear  theory,  similar  to  the  problem 
encountered  in  [86]. 

That  the  theory  works  very  well  in  the  cases  of  figures  1.4,  1.5,  1.6,  and  1.7 
but  not  for  figures  1.1  and  1.2  is  not  surprising-  when  the  stable  steady  state  is 
far  from  the  concentric  ring  solution,  the  linear  theory  is  less  likely  to  apply. 

Two-particle  minimizers  also  occasionally  break  symmetry  (in  the  sense  that 
the  steady  states  for  the  species  I  and  II  particles  differ  by  more  than  a  simple 
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reflection  or  rotation),  pictured  in  the  first  panel  of  figure  1.2.  Particle  densities 
from  figure  1.1  are  asymmetric  for  certain  parameters  (see  the  sixth  panel)  and 
may  be  supported  on  domains  with  irregular  geometry,  including  cusps.  Some 
simulations  show  dependence  on  initial  conditions,  which  did  not  manifest  for  the 
single  species  problem.  Of  course,  it  is  not  guaranteed  that  the  gradient  flow  (1.4), 
(1.5)  reach  a  global  minimizer  of  the  potential  (1.3).  This  minor  dependence  on 
initial  conditions  seems  to  indicate  that  the  energy  landscape  for  the  two-species 
problem  is  more  complex,  and  that  the  gradient  flow  occasionally  comes  to  rest 
at  local  minima. 

Simple  structures  directly  relevant  to  self-assembly,  such  as  alternating  particle 
rings  or  chains,  may  also  be  observed;  see  figures  1.4  and  1.8.  The  self-recognition 
and  separation  into  two  rings  replicates  the  phenomena  observed  in  [59],  and 
occurs  over  a  very  long  time  scale  (t  ~  1.4  x  105)  relative  to  that  of  the  formation 
of  the  rings  (t  ~  102). 

1.8  Conclusion 

Two-species  particle  aggregation  systems  have  a  rich  solution  structure,  including 
densities  with  rings,  spots,  and  radial  or  bilateral  TV-fold  symmetry  which  often 
concentrate  on  or  near  co-dimension  one  surfaces.  Considering  a  continuum  limit 
as  the  number  of  particles  tends  to  infinity  results  in  a  PDE  system  formulation 
similar  to  that  of  the  vortex  sheet  problem  [62,  78].  Linear  stability  analysis 
successfully  characterizes  the  steady  states  which  form,  verified  numerically  in 
section  1.7,  and  linear  well-posedness  of  the  PDE  system  is  considered  in  section 
1.6. 

Future  work  could  address  the  three  and  higher  dimensional  versions  of  the 
problem,  weakly  nonlinear  analysis  of  bifurcations  from  rings  to  other  steady 
states,  the  second-order  problem,  the  n-species  problem,  and  the  inverse  problem 
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lri  =  1  lr2  =  1  lr3  =  0.005 


Figure  1.8:  Alternating  particle  chains  arising  from  the  Morse  potential,  numbers 
of  particles  Ni  =  N2  =  N  with  N  —  8, 20,  80,  200, 400,  800  from  top  left  to  bottom 
right.  The  first  few  panels  show  that  the  particles  seem  to  form  effective  dipoles 
because  the  inter-species  repulsion  length  scale  is  so  small.  When  the  number 
of  particles  increases,  the  confining  nature  of  the  potentials  causes  them  to  pack 
closer  together  and  chains  form,  as  in  panel  5.  As  N  increases  further,  the  particles 
begin  to  form  a  two-dimensional  lattice  structure  (panel  6). 
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of  constructing  potentials  with  prescribed  instabilities  or  patterns  [86].  In  addi¬ 
tion,  the  two-species  system  allows  for  the  unique  possibility  of  nontrivial  H-stable 
ground  states  which  are  outside  the  scope  of  the  co-dimension  one  analysis  here; 
c.f.  [56]. 

Figure  1.9  shows  a  power  law  example  which  completely  leaves  a  co-dimension 
one  manifold,  and  seems  to  exhibit  an  effective  phase  separation  or  surface  tension 
arising  from  the  nonlocal  interactions.  As  the  inter-species  repulsion  singular¬ 
ity  becomes  weaker  than  the  intra-species  repulsion  singularity,  black  and  white 
particles  go  from  self-segregating  to  mixing.  When  the  inter-species  repulsion 
is  substantially  weaker,  exhibited  in  the  far  right  panel  of  figure  1.9,  a  regular 
alternating  lattice  structure  emerges  in  large  portions  of  the  steady  state.  The 
formation  of  lattices  in  the  single  species  problem  is  not  unfamiliar;  see  [39]  for 
similar  phenomena  in  the  single-species  problem.  The  first  two  panels  of  figure 
1.9  correspond  to  local,  not  global  minimizers  of  the  potential  1.3,  and  a  steady 
state  consisting  of  a  straight  line  interface  between  the  black  and  white  particles 
has  slightly  lower  energy.  The  steady  states  in  the  right  two  panels,  however, 
have  lower  energy  than  the  separated  black/white  half  disks.  The  analysis  carried 
out  in  this  chapter  applies  to  solutions  supported  along  one-dimensional  curves, 
and  as  such  does  not  apply  to  the  phenomena  in  figure  1.9;  however,  the  phase 
separation  effects  could  be  observed  on  a  co-dimension  one  surface  in  R3. 
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CHAPTER  2 


Sparse  Representations  for  Multiscale  PDE 


2.1  Preface 


This  chapter  discusses  a  new  area  of  research  (see  the  references  in  section  2.2) 
applying  sparse  modeling  techniques  to  the  numerical  solutions  of  PDE  with  mul¬ 
tiple,  separated  length  scales  [61].  A  typical  example  is 


due 


d 


dt  dx 


due 


dx 


— - — —  (  a(x/e)——  =  0  on  [0,  27t]  periodic, 


ue(x,  0)  =  u60(x),  a (x)  oscillatory 


where  e  is  near  zero.  The  highly  oscillatory  coefficient  a(x/e)  introduces  e-scale 
behavior  into  the  solution. 


The  reason  for  introducing  sparse  models  is  to  take  advantage  of  efficient  algo¬ 
rithms  available  for  sparse  data  structures.  The  results  presented  here  constitute 
advancements  in  numerical  analysis,  but  a  basic  background  in  sparse  modeling 
and  computation  with  sparse  data  will  help  put  them  in  context.  This  preface 
provides  that  background. 


2.1.1  Sparse  Representations  and  Sparse  Modeling 

A  vector  x  G  Mn  is  said  to  be  sparse  if  most  of  its  entries  are  zero.  A  signal  (or 
image  or  general  data  point)  b  G  Mm  is  said  to  admit  a  sparse  representation  with 
respect  to  a  possibly  overcomplete  basis  (or  dictionary )  D  G  Mmxn  if 

Da;  ~  b. 
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This  is  typically  meant  in  the  L 2  sense,  i.e.  ||Dx  —  b ||2  is  small.  A  sparse  model  of 
a  signal  or  set  of  signals  is  the  prior  knowledge  that  the  signals  admit  a  sparse  rep¬ 
resentation  with  respect  to  some  known  dictionary  D,  or  that  some  transformation 
of  the  signals  is  sparse.  Cases  in  which  the  signals  have  a  sparse  representation  are 
called  synthesis  models ;  those  in  which  a  transformation  of  the  signals  is  sparse 
are  called  analysis  models. 

Sparse  models  are  useful  because  they  provide  a  powerful  form  of  prior  knowl¬ 
edge  (this  statement  is  made  precise  in  Bayesian  terms  at  the  end  of  this  section) 
which  is  crucial  for  ill-posed  inverse  problems.  Sparse  models  are  particularly  pow¬ 
erful  because  they  are  widely  applicable  (many  types  of  data,  including  images 
and  signals,  have  sparse  representations)  and  because  they  can  be  used  to  formu¬ 
late  the  solutions  to  ill-posed  inverse  problems  involving  this  data  as  minimizers 
of  energy  functionals  which  are  computationally  tractable. 

To  illustrate  this  point,  consider  the  problem  of  removing  noise  from  a  signal 
which  is  assumed  to  have  some  unknown  sparse  representation  x  with  respect  to 
a  known  dictionary  D.  That  is,  the  true  signal  is  equal  to  D  t  for  some  sparse  x. 
If  the  noise  is  denoted  by  e,  then  we  observe  b  :=  Drr  +  e  and  wish  to  recover  Drc. 

The  sparse  model  is  what  makes  this  possible-  without  it,  the  problem  is 
completely  underdetermined.  We  expect  that  the  noise  is  not  too  large,  so  that 

||e|| 2  =  1 1 Da:  —  6|| \  is  small,  (2.1) 

and  also  that  x  is  sparse: 

||a;||0  is  small,  (2.2) 

where  ||a;||0,  the  L°  ‘norm’,  denotes  the  number  of  nonzero  entries  of  x.  (Note 
that  it  is  not  a  norm.) 

We  can  express  our  desire  to  find  a  representation  x  satisfying  both  require¬ 
ments  (2.1)  and  (2.2)  by  solving  an  optimization  problem  of  the  form 

x  =  argmin^.  ||Dx  —  6|| \  s.t.  ||a;||0  <  N 
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or 


x  =  argmin-r  ||a;||0  s.t.  ||Drc  —  6||  <  5. 

Alternatively,  the  search  for  x  can  be  cast  as  an  unconstrained  problem 

1  2 

x  =  arg  min  a,  A  ||x||0  +  -  ||Dx  —  b ||2  (2.3) 

where  A  which  is  a  parameter.  Solving  (2.3)  uses  the  L°  term  as  a  regularizer  to 
promote  sparsity  of  x,  and  balances  this  against  the  reconstruction  error  Drc  —  b. 

The  problem  with  (2.3)  is  that  the  L°  ‘norm’  is  totally  discontinuous  and 
so  the  resulting  optimization  problem  is  combinatorial  in  nature  and  NP-hard 
[64],  There  exist  efficient  greedy  algorithms  for  solving  (2.3)  with  performance 
guarantees  [84],  but  these  are  necessarily  approximate  in  nature.  Additionally,  it 
is  not  entirely  clear  that  the  L°  norm  is  the  best  way  of  measuring  sparsity:  its 
discontinuity,  at  least,  is  undesirable. 


An  alternative  to  approximately  solving  (2.3)  is  to  consider  its  convex  relax¬ 
ation 

1  2 

x  —  argmin-j,  A  +  -  ||Dx  —  &||2  .  (2.4) 

The  advantage  of  this  approach,  proposed  in  [21],  is  that  (2.4)  can  be  solved 
exactly  with  efficient  convex  optimization  methods  [40].  Another  possibility  is  to 
replace  the  L 1  norm  by  a  nonconvex  regularizer  [19,  20],  which  often  forfeits  the 
guarantee  of  finding  the  global  minimizer  but  performs  very  well  in  practice.  For 
the  purposes  of  this  introduction,  we  will  use  the  L 1  norm  with  the  understanding 
that  alternatives  are  available. 


This  approach  can  be  extended  to  handle  deconvolution,  upsampling  (or  su¬ 
perresolution),  and  other  inverse  problems  by  including  the  appropriate  forward 
operators  into  the  formulation  of  (2.4);  see  the  section  on  TV  image  reconstruc¬ 
tion  below.  All  that  is  required  is  that  a  reasonable  optimization  problem  can  be 

posed  with  balances  a  sparsity  term  (|| - 1| l5  say)  with  a  data  fidelity  term  (usually 

2 
2 
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Typical  examples  of  sparse  models  are  given  below: 


TV /ROF  Image  Model  [75]  This  is  the  only  analysis  model  discussed  here. 
The  prior  knowledge  is  that  images  (particularly  medical  images,  such  as  MRI  or 
X-ray  CT  scans)  are  approximately  piecewise  constant.  That  is,  their  gradients  are 
sparse.  An  image  x  can  thus  be  approximately  recovered  from  a  noisy  observation 
y  via 

1  2 

x  —  arg  min  x  A  ||  Va;^  +  -\\x  -  y\\2  . 

If  the  image  is  also  affected  by  an  operator  A  (which  could  correspond  to  blur, 
downsampling,  or  more),  the  optimization  problem  is 

1  2 

x  =  arg  min  x  A  ||  Vx\\1  +  -  \\Ax  -  y  ||2  . 

Transform-Domain  Models  [21]  and  Sparse  Coding  [67]  These  models 
assume  the  data  b  consists  of  an  image  or  image  patch  which  has  a  sparse  represen¬ 
tation  x  with  respect  to  either  a  wave  let /chirplet/warplet/shear  let/curve  let/etc 
dictionary  or  ‘learned’  dictionary  optimized  with  respect  to  a  given  set  of  signals. 
In  this  case,  x  represents  b  in  a  transform  domain  defined  by  the  dictionary  D. 
The  optimization  problem  generated  for  denoising  is  exactly  (2.4). 

Morphological  Component  Models  [77]  In  many  cases,  it  is  desirable  to 
separate  an  observed  signal  z  =  b\  +  b2  into  its  constituent  components  b\  and  b2. 
Examples  include  the  separation  of  speech  from  impulse  noise  and  cartoon/texture 
image  decomposition. 

Assuming  sparse  models  b\  «  and  62  ~  D2X2,  the  components  can  be 

separated  by  solving  the  optimization  problem 

„  „  „  „  1  „  ..  2 
xi,x2  =  arg  mm  Xl >X2  Ai  ||^i  ||1  +  A2  U^lli  +  -  +  D2x2  -  z ||2  • 
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Compressed  Sensing  [17,  18]  This  celebrated  field  combines  the  previously 
described  sparse  signal  models  and  efficient  convex  optimization  algorithms  with 
harmonic  analysis  and  random  matrix  theory  to  give  provable  guarantees  about 
the  recovery  of  transform-domain  sparse  signals  from  undersampled  measure¬ 
ments.  In  this  case,  the  transform  is  required  to  be  orthonormal  (such  as  a  number 
of  wavelet  and  discrete  Fourier  transforms). 

The  setup  is  as  follows:  we  suppose  that  /  €  Mn  is  a  signal  with  a  sparse 
representation  in  an  orthonormal  basis  \Er,  so  that  T'rr  =  /  with  x  sparse.  We 
observe  m  <  n  linear  functional  observations  of  the  form 

h  =  {(pi,  f)  =  (<Pi,  Vx) 

and  wish  to  recover  /.  The  Nyquist- Shannon  sampling  theorem  guarantees  that, 
in  general,  this  is  not  possible;  we  are  attempting  to  find  one  of  infinitely  many 
solutions  of  an  underdetermined  linear  system  A /  =  b.  However,  the  sparse  model 
T' x  =  f  provides  the  additional  information  necessary  to  reconstruct  /. 

In  other  words,  defining  A  as  the  matrix  with  (pj  ■  ■  ■  (p^  as  its  rows,  we  wish 
to  find  x  such  that 

A$a:  =  b  with  x  sparse. 

Along  the  lines  of  the  other  sparse  models  in  this  chapter,  we  propose  f  —  ^x 
where 

1 

x  =  arginine  A  +  -  ||A\Ihr  —  6||  . 

The  point  which  distinguishes  compressed  sensing  from  other  sparse  modeling 
contexts  is  that,  under  certain  conditions  on  A,  solving  the  above  optimization 
problem  is  guaranteed  to  succeed.  If  A  has  Gaussian  iid  entries  and  x  has  S 
nonzero  entries,  then  the  number  of  rows  (or  measurements)  required  is  just  m  = 
0{S  log  n)  [17]. 
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Parsimonious  Statistical  Models  [81]  and  Machine  Learning  L 1  regular¬ 
ization  appeared  in  [81]  as  a  method  for  enforcing  the  prior  knowledge  that  most 
of  the  coefficients  in  a  linear  regression  should  be  zero. 

The  standard  form  of  a  linear  regression  is  to  assume  that  an  output  y  sampled 
in  n  instances  y\i\,  i  =  1 . . .  n,  depends  on  inputs  X\[i\ . . .  xp[i]  in  a  form  which  is 
linear  in  the  parameters  /3\ ...  /3p: 

y[i]  =  f3\Xi[i]  +  . . .  +  f3pXp[i\  +  e[i\  for  all  i 

where  the  e[i]  ~  1V(0,  a2)  are  iid  zero- mean  Gaussian  random  variables.  In  matrix 
form,  the  model  can  be  written 

y^Xf3 

where  the  jth  column  of  X  is  Xj,  j  —  1 . .  .p.  Maximum  likelihood  estimation  for 
this  problem  leads  to  [43] 

13  =  argmm„  i  Hj,  -  X/J||  =  (XTXr'(XTt/). 

However,  if  n  <  p  the  situation  is  in  general  hopeless  because  the  data  y  could  be 
perfectly  explained  by  infinitely  many  parameter  choices  f3.  Even  if  n  >  p,  there 
is  a  danger  of  overfitting:  it  may  be  reasonable  to  expect  that  many  of  the  f3j  are 
in  fact  zero  because  the  corresponding  variables  Xj  have  only  a  negligible  effect 
on  y ,  but  the  maximum  likelihood  solution  will  have  fi.j  ^  0  and  use  these  Xj  to 
compensate  for  some  of  the  error  e. 

The  prior  knowledge  that  many  components  of  /3  are  zero  can  be  incorporated 
by  modifying  the  optimization  problem  so  that 

P  —  arg min  p  \  ^  \\y  -  X/9|| , 

known  as  LASSO  [81].  This  use  of  L 1  regularization  is  ubiquitous  in  modern 
data  mining,  especially  in  cases  where  there  are  thousands  of  input  variables  with 
negligible  effects.  The  LASSO  eliminates  them  automatically  with  the  tuning  of 
a  single  parameter  A. 
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In  [91],  the  authors  use  the  above  formulation  for  face  recognition.  The  x3  are 
example  (cropped  and  normalized)  faces  from  a  database,  and  y  is  a  face  to  be 
recognized.  The  coefficients  f3  are  then  either  thresholded  or  used  as  inputs  to  a 
support  vector  machine  or  other  classifier  to  determine  which  individual  in  the 
database  the  face  corresponds  to. 

Outlier  Models  and  Low-Rank  +  Sparse  Decompositions  [16]  In  many 
contexts,  such  as  collaborative  filtering  and  topic  modeling,  it  is  assumed  that  a 
matrix  Y  G  Mnxm  of  observations  can  be  approximately  factored 

Y  «  UV 

where  U  G  M.nxk  and  V  G  Mfcxm  with  k  much  smaller  than  either  m  or  n. 

However,  it  may  be  the  case  that  a  relatively  small  number  of  entries  of  Y 
do  not  fit  the  low-rank  model.  Along  the  lines  of  the  morphological  component 
model  discussed  previously,  we  can  formulate 

Y  «  UV  +  S 

where  the  majority  of  entries  of  S  are  zero.  This  leads  to  the  optimization  problem 

U,  V,  S  =  arg min U)V,s  i  ||UV  +  S  -  Y||J  . 

Other  regularizes,  such  as  the  nuclear  norm  [72],  can  be  used  to  limit  the  rank 
of  a  matrix  directly  rather  than  resorting  to  the  factorization  Y  «  UV.  The 
advantage  is  that  using  the  nuclear  norm  leads  to  a  convex  optimization  problem 
rather  than  the  nonconvex  one  above.  The  low-rank  %  sparse  decomposition  finds 
applications  in  video  foreground-background  separation  and  other  areas  [16]. 

Aside:  Sparsity  as  a  Bayesian  Prior 

In  the  above,  the  term  “prior”  has  been  used  loosely  to  describe  additional  knowl¬ 
edge  used  to  solve  ill-posed  inverse  problems,  separate  a  signal  into  components, 
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and  seek  meaningful  solutions  to  large  linear  regression  problems.  Here,  we  show 
that  this  usage  coincides  with  the  way  the  term  is  used  in  Bayesian  statistics. 

To  illustrate  this  point,  we  consider  the  problem  of  recovering  a  sparse  vector 
x  G  Mn  from  noisy  measurements  b  =  Ax  +  e,  where  A  G  Mnxm  and  e  is  a  vector 
of  mean-zero  Gaussian  noise  with  variance  a2 .  We  specify  a  Laplace  prior 

p(x)  =  C\  exp{—  Hxlh  /A} 

which  is  the  most  common  sparsity  prior  for  x,  but  there  are  other  choices.  C\  is 
a  normalizing  constant. 

The  likelihood  is  Gaussian: 

p{ Dx  —  b\x)  =  C2ex p  ||Dx  -  b\\l\  , 

and  so  the  posterior  is  given  by  Bayes’  rule  as 

f  1  2 

p(x\b)  oc  p(b\x)p(x)  oc  C\C2  exp  <  —  /A  —  — ^  | Da:  —  b\\2 

The  posterior  attains  its  maximum  at 

x*  =  argmin-r  —  H —  ||Dx  —  6|| \ 

A  2 

which  is  exactly  the  optimization  problem  proposed  above.  In  Bayesian  terms, 
using  this  value  of  x  at  the  mode  of  the  posterior  is  called  maximum  a-posteriori 
(MAP)  estimation.  A  true  Bayesian  treatment  would  avoid  using  a  point  estimate 
for  x  at  all,  preferring  the  whole  posterior  distribution  instead.  When  a  point 
estimate  is  necessary,  the  usual  Bayesian  choice  is  the  mean  of  the  posterior  instead 
of  its  maximum.  Computing  the  posterior  mean  is  usually  done  with  Markov 
Chain  Monte  Carlo  (MCMC)  methods,  or  variational  inference. 


2.1.2  Computation  with  Sparse  Data 

The  application  to  numerical  PDE  we  are  considering  does  not  involve  any  ill- 
posed  problems  for  which  sparsity  provides  the  necessary  additional  information 
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Figure  2.1:  Sparsity  pattern  of  a  typical  finite  element  matrix.  Only  the  black 
pixels  correspond  to  nonzero  entries.  Figure  courtesy  of  Wikipedia  [88]. 

to  recover  a  solution.  The  PDEs  we  consider  are  well-posed,  and  have  unique 
solutions. 

Instead,  the  primary  role  of  sparsity  in  this  chapter  is  to  allow  for  efficient 
computation.  This  is  already  common  practice  in  numerical  analysis  and  scientific 
computing,  and  it  is  one  of  the  primary  advantages  of  iterative  methods  such  as 
Conjugate  Gradient,  MINRES,  and  GMRES  [41]. 

The  basic  idea  is  that,  when  most  of  the  entries  in  a  matrix  or  vector  are 
zero,  memory  and  time  can  be  saved  by  storing  only  the  location  and  value  of  the 
nonzero  entries.  For  example,  matrices  arising  in  finite  element  methods  have  the 
property  that  the  number  of  nonzero  entries  in  each  row  (and  column)  is  bounded 
by  the  maximum  number  of  elements  adjacent  to  any  given  element,  which  stays 
fixed  as  the  mesh  is  refined.  Figure  2.1  shows  such  a  matrix. 

If  a  sparse  matrix  A  is  size  N  x  N,  naive  multiplication  against  an  ordinary 
dense  vector  requires  0(1V3)  floating  point  operations.  However,  because  A  is 
sparse,  there  exist,  say,  only  k  nonzeros  per  row.  By  only  multiplying  the  nonzero 
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entries  of  A  times  a  vector,  the  number  of  floating  point  operations  can  be  reduced 
to  0(kN).  Since  k  stays  fixed  as  N  increases,  this  is  a  huge  complexity  gain. 


2. 1.2.1  Efficient  Convolution  of  Sparse  Vectors 


This  chapter  focuses  on  solving  PDEs  of  the  form 

due  d  (  due\ 

&)=°  °M0, 2*]  periodic, 

ue(x,  0)  —Uq(x),  f(x/e)  oscillatory 

in  the  Fourier  domain.  Explicit  methods  and  implicit  methods  based  on  Krylov 
subspace  solvers  all  require  repeatedly  applying  the  elliptic  operator 

Yx(fix/e)%)' 


which  comes  down  to  convolutions  of  the  form 


x*y, 

where  x  and  y  are  sparse,  length- N  vectors  representing  /  and  u  discretized. 

If  x  has  nx  nonzero  entries  and  y  has  ny ,  then 

N—l 

(x  *  y)\i\  =  ^  x[i  -  j  mod  N]  ■  y\j] 

3=0 

can  be  thought  of  as  the  summation  of  ny  scaled  shifts  of  x ,  so  that  computing 
the  convolution  is  tantamount  to  adding  together  ny  sparse  vectors  (it  can  also 
be  thought  of  as  a  matrix- vector  multiplication).  We  assume  that  x  and  y  are 
stored  as  lists  of  (index,  value)  pairs,  the  most  commonly  used  format  for  sparse 
vectors.  Specifically,  we  assume  there  is  a  queue  data  structure  qx  which  represents 
x.  That  is,  qx  has  nx  entries  {qx[k]}^= o\  each  of  which  is  an  {index,  value)  pair 
{qx[k}. first,  qx[k\. second)  such  that  x[qx[k}. first]  =  qx[k\. second.  We  assume  that 
the  queue  allows  iteration  over  its  elements,  and  0(1)  insertion  at  the  back,  and 
0(1)  deletion  from  the  front  (these  are  all  standard  assumptions  for  queues).  For 
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example,  such  a  queue  could  be  implemented  easily  with  a  linked  list.  We  assume 
qy  is  the  analogous  queue  representing  y. 

For  example,  if  x  is  a  vector  of  length  10  with  three  nonzero  entries  a;[l]  =  — 1, 
x[5]  =  1,  and  x\7]  =  6,  the  queue  qx  would  consist  of  the  pairs  (1,  —1),  (5, 1),  (7,  6). 


We  now  discuss  several  possible  algorithms  for  computing  the  convolution  of 
two  sparse  vectors,  which  is  the  main  bottleneck  in  our  numerical  procedure. 

Merge  Approach  If  we  assume  that  the  queues  qx  and  qy  representing  x  and 
y  are  sorted  by  increasing  index  (as  in  the  above  example),  the  problem  can  be 
viewed  as  one  of  merging  sorted  lists  into  a  longer  sorted  list.  When  there  is  a 
‘tie’,  i.e.  two  of  the  vectors  to  be  added  have  a  nonzero  in  the  same  index,  the 
corresponding  values  are  simply  added  in  the  result. 

Naively  merging  the  sorted  lists  by  searching  for  the  minimum  index  among 
the  ny  shifts  of  x  would  result  in  an  algorithmic  complexity  of  0(riynx).  At  each 
iteration,  the  front  entry  of  each  of  the  ny  shifts  would  be  checked  to  determine 
the  minimum,  a  0(ny)  cost.  Each  iteration  adds  only  one  additional  coefficient 
to  the  result,  so  nxny  total  iterations  are  required,  leading  to  the  overall  0[;ri^nx) 
complexity.  This  has  the  advantage  of  being  independent  of  the  total  grid  size  N, 
but  as  it  turns  out  this  algorithm  is  far  from  optimal. 

Keeping  track  of  the  minimum  index  during  the  merge  can  be  accomplished 
more  efficiently  with  a  heap  data  structure  [89],  which  is  ideally  suited  for  this 
purpose.  Details  of  the  algorithm  revolve  around  data  structures  rather  than 
mathematics,  and  a  more  efficient  algorithm  will  be  presented  next,  so  pseudocode 
for  this  approach  is  omitted.  The  complexity  is  0(nxny  log (ny)),  a  substantial 
improvement  over  the  naive  algorithm. 
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Index-tracking  Approach  Ideally,  it  would  be  possible  to  compute  the  con¬ 
volution  (2. 1.2.1)  with  0(nxny)  complexity.  This  is  the  best  possible,  because 
the  sum  in  (2. 1.2.1)  consists  of  nxny  terms  which  each  require  a  floating  point 
multiplication. 

It  turns  out  that  this  lower  bound  is  achievable,  and  the  key  to  reaching  it  is 
to  notice  that  the  extra  log(ny)  factor  in  the  complexity  of  the  merge  approach  is 
due  to  the  computational  cost  of  keeping  track  of  the  minimum  index.  However, 
this  can  be  avoided  in  two  ways.  The  first  is  to  use  the  fact  that  the  vectors  to 
be  added  are  not  arbitrary,  but  rather  shifts  of  a  single  vector  x.  Therefore,  the 
vectors  which  will  contribute  to  a  nonzero  entry  in  a  particular  index  of  the  result 
can  be  known  ahead  of  time  rather  than  computed  with  the  heap. 

The  second  approach,  which  is  more  general,  is  to  use  an  auxiliary  vector 
to  store  the  result  with  an  additional  array  to  keep  track  of  which  entries  are 
nonzero.  This  way,  the  indices  of  the  nonzero  entries  do  not  need  to  stay  sorted. 
Our  pseudocode  for  the  algorithm  assumes  the  following  notation: 

•  qx  is  a  queue  representing  x  in  (index,  value)  pairs  as  described  above. 

•  qy  is  an  analogous  queue  representing  y. 

•  z  is  an  empty  queue  of  the  same  form,  which  will  store  the  result  of  the 
convolution. 

•  ind_list  is  an  empty  queue,  which  will  hold  integers  in  the  range  0  :  N  —  1. 
Let  ind_list. length  denote  its  length. 

•  t  is  an  array  with  entries  Recall  that  arrays  offer  0(1)  entry 

lookup,  which  will  be  crucial  later. 

Pseudocode  for  the  sparse  convolution  -  index-tracking  algorithm  is  given  be¬ 
low,  and  accomplishes  MAX  ITERATIONS  sparse  convolutions  with  shrink. 
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Sparse  Convolution  -  Index-tracking  Algorithm 

initialize  t[i]  =  0  for  all  i  =  0  :  N  —  1  >  O(N)  cost 

for  counter  =  0:M AX  ITERATIONS  —  1  do 
for  j  —  0 :  ny  *-  1  do 
for  k  —  0 :  nx  —  1  do 

if  (t[qx[k\. first  +  qy[j]. first  mod  N]  =  0)  then 

append  qx[k\. first  +  qy[j]. first  mod  N  to  ind_list  >  0(1)  cost 

end  if 

t [qx[k\- first *jr  qy[j], first  mod  N]  +=  qx[k\. second  x  qy[j]. second  > 

0(1)  cost 

end  for 
end  for 

for  i  =  0  \ind_list. length  —  1  do 
if  (\t[ind_list[i]\\  >  A)  then 

append  (ind_list[i\,shrmk(t[ind_list[i]]))  to  z  >  0(1)  cost 

end  if 

set  t[ind_list[i]]  =0  >  0(1)  cost 

end  for 
end  for 
return  z 
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The  nested  loops  over  j  and  k  lead  to  a  total  cost  of  0(nxny),  the  best  one 
could  hope  for.  At  each  step  of  this  loop,  at  most  one  entry  is  added  to  ind_list, 
so  the  loop  over  %  costs  at  most  0(nxny )  as  well.  Thus  the  whole  cost  is 

O  (nxnyM  AX —ITERATIONS)  +  O(N) 
where  the  only  0(N )  cost  is  the  initialization  in  the  first  line,  which  is  cheap. 

This  algorithm  strengthens  the  results  of  the  chapter  by  demonstrating  that 
efficient  computation,  independent  of  the  grid  size  N  after  the  first  iteration,  is 
possible.  There  are  other  approaches  to  multiscale  PDE  with  scale  separation  (e.g. 
[28])  which  leverage  sparse  representations  of  the  solution  for  efficiency,  but  so  far 
none  with  algorithmic  complexity  totally  independent  of  the  grid  size  (except  for 
the  initial  and  final  Fast  Fourier  Transforms). 

2.2  Introduction 

Partial  differential  equations  with  multiple  length  scales  are  fundamental  to  mod¬ 
eling  various  physical  problems  including  composite  materials,  wave  propagation 
in  inhomogeneous  media,  crystalline  solids,  and  flows  with  high  Reynolds  number 
(fluid  mechanics).  Typically,  these  problems  involve  a  wide  range  of  scales,  with 
each  scale  corresponding  to  a  level  of  physical  processes.  However,  in  some  cases, 
the  problem  is  scale  separable,  in  the  sense  that  the  mathematical  representation 
of  the  dynamics  involve  one  fine  scale  and  one  course  scale.  Even  in  this  case, 
accurate  numerical  methods  for  solving  these  PDE  can  be  computationally  ex¬ 
pensive  since  resolving  both  the  coarse  and  fine  scales  simultaneously  requires  a 
spatial  resolution  dominated  by  the  fine  scale. 

Over  the  past  decades,  various  approaches  have  been  taken  to  overcome  this 
difficulty.  In  some  cases,  it  is  possible  to  derive  an  asymptotic  approximation  for 
the  effect  of  small  scales  on  the  solution  [70].  When  this  is  not  possible,  many 
other  techniques  have  been  proposed.  A  multiscale  finite  elements  method  can 
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be  used  to  solve  linear  elliptic  homogenization  equations  (see  [47]),  and  has  found 
many  applications  to  other  multiscale  problems.  The  equation-free  methods  use 
accurate  small  scale  and  short  time  solvers  to  capture  fine  scale  behavior  and 
use  them  to  govern  the  related  course  scale  behavior  [49].  The  heterogeneous 
multiscale  method  [34]  is  a  general  numerical  approache  which  also  uses  the  scale 
separation  of  the  problem  to  generate  solvers  on  the  micro  and  macroscopic  levels. 
In  [66],  a  projection  based  approach  is  used  to  construct  an  adaptive  multiscale 
algorithm  for  elliptic  homogenization  equations.  And  more  recently,  a  sparse 
transform  method  [28]  exploits  the  scale  separability  of  linear  homogenization 
problems  to  construct  a  fast  direct  solver.  The  body  of  literature  on  multiscale 
models  is  large,  and  we  only  mention  some  of  the  popular  methods.  For  more 
detail  on  general  numerical  methods  for  multiscale  problems  see  [35,  34]  and  the 
citations  therein. 

In  this  work,  we  will  focus  our  attention  on  linear  partial  differential  equations 
with  multiscale  behavior  either  in  the  medium  or  in  a  source  term.  Following  the 
work  of  [76],  which  used  an  L 1  optimization  method  to  compress  the  Fourier  coef¬ 
ficients  of  the  solution,  we  build  efficient  solvers  for  periodic  multiscale  problems. 
In  particular,  we  will  use  the  sparse  Fourier  structure  of  solutions  to  construct 
numerical  methods  which  solve  the  problem  directly,  without  the  separating  the 
micro  and  macro  scales  explicitly. 

L 1  optimization  and  its  related  models  are  at  the  center  of  many  problems  in 
the  fields  of  imaging  science  and  data  analysis,  see  for  example  [17,  16,  31,  15]. 
Due  to  the  connection  with  sparse  models  for  compressive  sensing,  recent  works 
have  introduced  L 1  techniques  for  numerical  partial  differential  equations.  For 
example,  in  [76]  L 1  regularized  least  squares  was  used  to  sparsely  approximate  the 
Fourier  coefficients  in  multiscale  dynamic  PDE  (and  in  this  work  we  expand  that 
approach).  In  [65,  68,  69],  eigenfunctions  with  compact  support  were  constructed 
to  efficiently  solve  problems  in  quantum  mechanics.  Also,  in  [46]  an  L 1  nonlinear 
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least  squares  model  was  used  to  sparsely  recover  coefficients  of  a  second  order 
ODE  which  are  related  to  constructing  intrinsic  mode  functions.  In  [11],  low- 
rank  libraries  are  used  to  sparsely  approximate  solutions  to  dynamical  systems 
and  thereby  identify  bifurcation  regimes.  Some  theoretical  results  are  provided  in 
[14]  for  PDE  with  L1 -terms,  related  to  some  of  these  models.  For  more  detailed 
analytic  results,  see  [8,  9,  10]  which  laid  the  theoretical  groundwork  for  these 
equations. 

In  this  chapter,  we  continue  the  work  of  [76]  to  leverage  the  sparsity  of  solutions 
in  order  to  design  an  efficient  numerical  scheme.  However,  we  also  impose  sparsity 
of  the  update  operator  to  improve  the  complexity  while  retaining  a  similar  level 
of  accuracy.  We  show  some  theoretical  results  for  the  sparse  spectral  scheme  and 
sparse  operator-sparse  solution  spectral  scheme.  In  particular,  we  provide  error 
bounds  between  the  solution  and  the  sparse  approximation  as  well  as  complexity 
bounds  on  the  algorithm.  Also,  we  continue  to  make  connections  between  L 1 
based  methods  and  multiscale  problems  through  a  denoising  interpretation  of  the 
homogenization  expansion  of  the  solution. 

The  outline  of  this  work  is  as  follows.  In  Section  2.3,  we  recall  the  explicit 
scheme  from  [76]  and  in  2.4  propose  an  implicit  version  as  well  as  a  sparse  operator 
approximation.  Theoretical  results  are  provided  in  Section  2.5.  A  discussion  on 
well-posedness  is  given  in  Section  2.6  and  a  denoising  interpretation  of  the  method 
is  given  in  Section  2.7.  In  Section  2.8,  some  algorithmic  analysis  is  provided.  The 
algorithm  is  tested  on  numerical  examples  in  Section  2.9,  with  concluding  remarks 
given  in  Section  2.10. 

2.2.1  Notation 

•  a  -  (or  A  for  anisotropic  problems)  the  medium  inhomogeneity,  a!  is  the 
sparse  approximation  of  a. 
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•  n  -  the  shrink  size  variable.  //'  is  the  corresponding  variable  for  sparse 
operator  approximation. 

•  k  -  the  Fourier  space  variable,  with  positive  and  negative  frequencies. 

•  Q  -  either  a  general  numerical  scheme  or  the  matrix  corresponding  to  a 
one-step  linear  numerical  scheme. 

•  L  an  elliptic  operator.  L  is  the  operator  when  applied  in  the  Fourier 
domain  and  L ^  is  its  discretization.  L'h  is  the  sparse  approximation. 

2.3  Preliminary 

We  will  consider  linear  multiscale  problems  where  the  solutions  are  sparse  in  the 

Fourier  domain  [28,  76].  For  example,  consider  the  parabolic  problem: 


due 

~dt 


d_ 

dx 


0  on  [0,  27r]  periodic 


ue(x,  0)  =  ue0(x),  a(x/e)  oscillatory. 


(2.5) 


Figure  2.2  shows  the  solution  in  physical  and  Fourier  space.  This  phenomenon  is 
common  in  multiscale  PDE:  distinct  length  scales  manifest  strikingly  as  sparsity 
in  the  frequency  domain. 

To  compute  solutions  which  are  truly  sparse  in  the  frequency  domain  (and  not 
just  approximately  sparse  with  many  noisy  small  magnitude  coefficients),  it  was 
proposed  in  [76]  to  solve  an  7 '-regularized  least  squares  problem  to  obtain  a  sparse 
approximation  of  u  (the  Fourier  transform  of  u ).  We  summarize  the  method  here. 

Given  numerical  iterates  un, . . .  ,un~q  and  a  numerical  update  scheme  of  the 
form 

un+x  =  Q(un,...,un~q), 
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Figure  2.2:  Left:  Solution  of  (2.5)  with  Fourier-sparse  initial  data  in  physical 
space.  The  small  rectangle  shows  the  axis  limits  of  the  zoomed  in  plot  to  the 
right.  Right:  Zoomed  in,  showing  fine  scale  oscillations.  Bottom:  solutions  in 
Fourier  space  (the  y-axis  for  all  Fourier  space  plots  is  on  a  log10  scale).  Of  the 
N  =  2048  Fourier  coefficients,  only  153  have  magnitude  larger  than  10-10. 
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the  scheme  is  modified  by  defining  the  auxiliary  variable  v  =  Q[un, . . .  ,un  q )  and 
solving 


un+1  =  argminw  n  ||io|| x  +  -  || w  —  v\\\  , 


(2,6) 


where  the  d1  norm  for  complex  arguments  w  is  =  )>T  \wi\,  where  |  •  |  denotes 
magnitude.  Note  that  the  d1  norm  is  taken  in  the  Fourier  domain  and  not  physical 
space. 


For  a  one-step  linear  updating  scheme,  equation  (2.6)  can  be  written  as 


u 


71+1 


arg  ruin  w 


+  \\\w~  Qun 


2 

2 


where  Q  is  the  matrix  which  advances  the  discretized  solution  forward  in  time.  L 1 
regularized  least  squares  is  amenable  to  a  number  of  efficient  solution  methods, 
e.g.  [40].  The  problem  can  also  be  generalized  to  any  basis  or  overcomplete 
dictionary,  but  we  restrict  our  attention  to  Fourier  modes.  In  fact,  due  to  the 
orthogonality  of  the  Fourier  modes,  equation  (2.6)  decouples  and  the  minimizer 
can  be  given  exactly: 


1) 

un+1  =  shrink (v,  fi)  :=  max(|h|  —  /i,  0)-^-. 

M 

For  a  concrete  example,  the  forward  Euler  method  applied  to  =  Jk  [ a(x / 
has  the  form: 


un+l  =un  +  dtika*  (ik un ) 


where  k  is  the  wave  number  and  *  represents  convolution.  This  becomes: 
un+1  =  shrink  (un  +  dtika*  (i  k  un) ,  /i ) 

in  the  sparse  spectral  form.  By  exploiting  sparsity  in  the  frequency  domain,  the 
proposed  method  can  rely  on  sparse  data  structures  to  allow  for  high  resolution 
with  faster  numerical  simulations. 
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2.4  Proposed  Methods 


In  this  section,  we  will  discuss  two  new  extensions  of  the  sparse  spectral  method, 
namely  an  implicit  version  and  a  sparse  operator/sparse  solution  version.  Each 
come  with  their  own  advantages,  which  we  will  analyze  in  subsequent  sections. 

2.4.1  Implicit  Variation 

For  many  classes  of  problems  at  high  spatial  resolution,  explicit  schemes  are  im¬ 
practical  due  to  the  severe  time  step  restriction  required  for  stability.  We  can 
construct  an  implicit  scheme  for  the  problems  we  are  considering,  which  avoids 
these  restrictions  at  the  expense  solving  a  more  complex  L 1  problem  at  each  time 
step. 

Consider  the  general  linear  parabolic  equation  ut  +  Lu  =  f  with  schemes  of 
the  form 

Qun+ 1  =  un  +  dtfh. 

The  simplest  implicit  version  is  backward  Euler: 

(/  +  dtLh)un+1  =  un  +  dtfh ■, 

where  L  denotes  the  representation  of  L  in  the  Fourier  basis,  Lh  denotes  the 
discretized  version  of  this  operator  with  respect  to  a  grid  size  h  >  0,  and  fh 
denotes  the  Fourier  transform  of  /  sampled  at  the  corresponding  grid  points. 

For  a  scheme  of  this  form,  the  analogue  of  equation  (2.6)  is 

1  -  2 

un+l  —  argminw  ^11^1^  +  -  Qw  —  (un  +  dtfh)  (2.7) 

which  does  not  have  a  simple  explicit  representation.  In  addition,  the  optimality 
condition  for  Equation  (2.7)  requires  inverting  the  matrix  QTQ  which  will  often 
be  badly  conditioned. 
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When  L  is  a  uniformly  elliptic  operator,  the  eigenvalues  of  Q  —  I  +  dtLh  are 


positive  and  so  we  can  instead  consider  the  sparse  scheme  defined  by 

1 


/v  T7  -I- 1  • 

u  —  arg  mm  w 


H  IMIi  +  2W  Qw  ~  w  n  +  dtfh). 


(2.8) 


Similarly  for  time- independent  problems,  i.e.  Lu  =  /,  the  corresponding 
energy  is 

u  =  arg  min  w  /i  ||tt;||  x  +  -wT Lh.w  —  vo1 fh- 

Note  that  when  /x  —  0,  this  is  the  standard  variational  principle  for  elliptic  oper¬ 
ators.  We  will  see  that  solving  the  implicit  schemes  with  the  L 1  term  directly  is 
often  too  slow  to  be  practical.  The  reason  is  that  directly  applying  this  variational 
principle  to  find  the  solution  does  not  use  the  fact  that  the  solution  is  sparse  in 
order  to  speed  up  computations.  However,  in  Section  2.8.1  we  will  show  that  it  is 
possible  to  construct  an  efficient  algorithm  for  approximately  solving  the  resulting 
optimality  condition  arising  from  equation  (2.8). 


2.4.2  Sparse  Operator  Approximation 

For  uniformly  elliptic  linear  operators,  for  example  of  the  form 

Lu  =  —  div(a(x,  i/e)V«), 
the  standard  spectral  discretization 

Lhu  —  ka*  (ku) 

requires  a  convolution  at  each  iteration,  which  can  be  costly  even  when  u  is  sparse. 
However,  because  the  diffusion  coefficient  a  is  scale  separated,  we  can  define  a 
sparse  approximation  of  Lh  by 

L'hu  —  ka!  *  (ku) 

where  a!  is  a  sparse  approximation  of  a.  We  can  choose  a!  to  solve 

a  =  arg  min  w  /a'  ||io|| x  +  -  ||w  —  a\\l 
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which  again  results  in  a  closed  form  solution  given  by  the  soft  thresholding  a'  = 
shrink(a,  //).  An  alternative  is 

a  =  argminw  p!  ||w||0  +  -  ||w  —  a|| \ 

where  the  L°  ‘norm’  ||-||0  counts  the  number  of  nonzero  entries.  In  this  case, 
the  solution  is  given  by  hard  thresholding  a —  setting  all  coefficients  smaller  in 
magnitude  than  some  threshold  equal  to  zero. 

Soft  thresholding  is  contractive  and  benefits  from  many  desirable  smoothing 
properties  which  make  it  preferable  for  the  sparse  approximation  of  the  solution, 
which  will  be  discussed  below  in  Section  2.7.  For  a  sparse  approximation  of  the 
operator,  the  benefits  of  a  particular  choice  of  thresholding  are  less  clear  and 
therefore  we  consider  both. 

2.5  Theoretical  Remarks 

The  compressive  spectral  method,  or  sparse  scheme,  inherits  many  appealing 
properties  of  the  underlying  numerical  method  it  approximates.  In  general,  it  is 
at  least  as  stable  as  the  original  scheme  and  retains  the  order  of  accuracy. 

2.5.1  Contraction  and  Linear  Convergence 

The  following  two  theorems  show  that  the  explicit  and  implicit  numerical  schemes 
are  contractive.  This  result  is  similar  to  those  found  in  [9]. 

Theorem  2.5.1.  For  the  explicit  scheme  generating  time  steps  by 

<  un+1  =  shrink((I  —  dtLh)un  +  dtf,  p), 

if  ||J  —  dtLh\\op  <  1  then  the  iterations  are  contractive:  i.e.  ||un+1  —  un ||2  < 

||  Un  —  Un~  1  1 1 2  - 

Here  ||-||op  denotes  the  £2  operator  norm,  or  largest  singular  value. 
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Theorem  2.5.2.  For  the  implicit  scheme,  if  Lh  is  positive  semidefinite,  then  the 
iterations  are  contractive,  ||-uri+1  —  un H2  <  || un  —  un~l\\-2,  for  all  dt  >  0. 

The  proofs  of  these  two  theorems  are  similar,  and  reside  in  the  appendix. 

The  method  is  also  convergent.  In  particular,  for  the  correct  scaling  of  p,  we 
have  the  following  theorem. 

Theorem  2.5.3  (Linear  Convergence,  Explicit  Scheme).  Let  S  denote  a  linear 
spectral  numerical  update  scheme,  generating  time  steps  as 

un+1  =  Q(un,...,un~k), 

and  let  5)t  denote  the  spectrally  sparse  scheme,  which  generates  time  steps  as 

un+l  =  shrink{Q{ul , . . . ,  un~k),tf). 

Then  if  S  is  consistent  and  stable  (and  hence  converges),  and  if  p  =  0(dt1+s )  for 
some  5  >  0,  then  the  compressive  scheme  S 'M  converges.  If  p  =  0(dtp )  with  p  at 
least  the  order  of  the  local  truncation  error  of  S,  then  the  order  of  convergence  of 
S  is  not  impacted. 

For  the  implicit  scheme,  the  analogous  theorem  is  the  following. 

Theorem  2.5.4  (Linear  Convergence,  Implicit  Scheme).  Let  S  denote  an  implicit 
linear  spectral  numerical  update  scheme  for  the  PDE  ut  +  Lu  =  /  on  a  domain 
discretized  with  N  grid  points,  generating  time  steps  as 

(/  +  dtLh)un+1  =  un  +  dtfh , 

and  let  St,  denote  the  spectrally  sparse  scheme,  which  generates  time  steps  accord¬ 
ing  to  (2.8): 

1 

un+1  =  argminw  p  ||xu|| x  +  -U!T(/  +  dtLh)w  —  tcT(n))  +  dtfh). 
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Then  if  S  is  consistent  and  stable  (and  hence  converges),  and  if  n  =  0(dt1+s )  for 
some  5  >  0,  then  the  spectrally  sparse  scheme  Stl  converges.  If  p  =  0(dtp )  with  p 
at  least  the  order  of  the  local  truncation  error  of  S,  then  the  order  of  convergence 
of  S  in  L 2  is  not  altered. 

The  proofs  can  be  found  in  the  Appendix. 

2.5.2  Sparse  Operator  Approximation:  Implicit  Solver 

We  now  consider  the  error  incurred  by  the  sparse  operator  approximation  proposed 
in  Section  2.4.2.  The  continuum  case  is  discussed  in  detail,  and  the  proof  for  the 
case  of  discretized  operators  is  completely  analogous. 

The  usual  discretization  in  Fourier  space  of  a  general,  anisotropic,  divergence 
form  elliptic  operator 

Lu  =  — div(A(a;)Vu)  +  b(x )  •  V«  +  c(x)u 

results  in  a  matrix  (corresponding  to  convolution)  which  is  dense.  However,  it  is 
still  approximately  sparse  when  the  coefficients  A  and  b  are.  Approximating  A 
and  b  by  A'  and  b'  which  are  truly  sparse  in  Fourier  space  yields  an  operator  which 
is  far  more  efficient  to  store  and  to  work  with,  but  incurs  some  error.  Theorem 
2.5.5  quantifies  this  error. 

Theorem  2.5.5.  Let  u\  and  U2  be  solutions  to 


— div(AiVwi)  +  bi  ■  V«i  +  C\U\  =  /  (2.9a) 

-div(A2Vu2)  +  h  ■  Vit2  +  c2u2  =  f.  (2.9b) 


on  a  domain  Oc  with  periodic  boundary  conditions  and  the  constraint 


u2  =  0. 


'CL  JCL 


u  i  = 
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Require  also  that 


wTAiW  >  A  ||ic||2 , 

Q  -  ^div(6i)  >  0 

for  i  =  1,2.  Then 

\\ui  -u2\\Hi  <  UA_2(  d max  (Ai)^  -  ( A2)ij  + 

V  ^  1 

C  max  (bi)i-(b2)i  ^  +  C2  j|ci  —  C2||x  'j  ||/||2 

where  C  =  C(Q)  is  the  constant  from.  Poincare’s  inequality  [38],  and  the  Fourier 
series  of  the  matrices  and  vectors  Ai  and  bi  are  taken  entry-wise. 

This  form,  in  terms  of  (Ai)^  —  {A2)ij  ,  (&i);  —  (b2)i  ,  and  ||c!  —  is 

particularly  useful  because  the  coefficients  will  be  approximated  in  Fourier  space. 
The  reader  familiar  with  energy  estimates  for  elliptic  equations  will  see  that  the 
requirements  of  the  theorem  are  not  the  most  general  possible,  and  the  proof 
can  be  modified  to  handle  other  cases  individually  when  different  estimates  are 
desired. 

Proof.  Subtracting  the  first  equation  of  (2.9)  from  the  second,  then  adding  and 
subtracting  A{\/u2,  b\Wu2,  and  C\U2  gives 

-div(Ai Vic)  -  div[(Ai  -  A2)  Vu2]  +  h  ■  Vic  +  (bi  -  b2)Vu2  +  cxw  +  (cq  -  c2)u2  =  0, 
and  after  multiplying  by  w  and  integrating  by  parts,  we  glean 

A  J  |  Vw|2  dx  <  J  VicT A\Vw  +  ^ci  —  ^div(6i)^  w2  dx 

<  || Ali  -  A2\\op  1 1 V 1 1 2  ||VH|2  +  ll&i  -  ML  HVttalla  IMI2  +  ••• 

Hci-calLWalkL- 

Using  Poincare’s  inequality  and  ||^4i  —  A2\\op  <  d\\A±  —  v42 1| ^ , 

A  /  |  Vic  | 2  dx  <{d\\A1^A2\\0B  +  C  \\h  -  ML  +  U2  ||Cl  -  ML)  ||  Vu2||2 1|  Vtc||2 

Jn 
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and  thus 


l|V<0||2  <  A -‘(d  \\Ai  -  ^IL  +  C  lift!  -  ft2|L  +  C2  ||ci  -  e2||J  ||V«2||2 .  (2.10) 

Similarly  multiplying  the  equation  for  u2  by  u2 ,  integrating  by  parts,  and  applying 
Poincare’s  inequality  yields 

||Vu2||2<CA-1||/||2. 

Substituting  this  into  (2.10)  and  using  Poincare’s  Inequality  again,  we  get 

IK  -  yy  <  c\-\d ||A!  -  a2\\x + c lift!  -  ft2y  +  c2 iid  -  c2iy  h/ii2 . 

The  form  stated  in  the  theorem  follows  after 

Plloo  =  max  Pulloo  ^  max  Aij 

1,3  1,3  1 

and  the  analogous  inequality  with  b.  □ 

In  practice  memory  is  not  a  concern  due  to  the  convolutional  structure  of 
the  matrix  Lh  representing  an  elliptic  operator  in  Fourier  space,  but  the  sparse 
structure  of  the  operator  dramatically  reduces  computation  complexity  (Section 
2.8.2). 

2.5.3  Sparse  Operator  Approximation:  Explicit  Solver 

The  discrete  analogue  of  Theorem  2.5.5  covers  numerical  schemes  with  implicit 
time  steps,  each  of  which  require  solving  an  elliptic  problem  with  a  sparsely  ap¬ 
proximated  operator.  Effectively,  it  allows  us  to  estimate 

IP”1  -  F~'IU 

where  P  is  a  sparse  matrix  approximating  that  of  the  discretized  full  elliptic 
operator  Q ,  and  ||-||  refers  to  the  L 2  matrix  operator  norm,  or  largest  singular 
value.  On  the  other  hand,  for  explicit  schemes,  we  are  concerned  about 

\\Q-P\\op 
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which  we  will  consider  directly. 


Theorem  2.5.6.  Let  L  be  the  elliptic  operator  defined 

Lv  =  —  div(aVc) 

and  let  Q  be  its  Fourier  discretization 

Qu  —  kdh*  ( ku ) 

where  k  denotes  the  vector  of  Fourier  mode  frequencies  and  ah  is  the  discretized 
domain  inhomogeneity  coefficient  in  the  elliptic  operator.  Analogously,  let 

Pu  =  k  a'h  *  (ku). 

Then 

\\Q  -  P\\o„  <  K2  \\*k  - 

where  K  is  the  highest  frequency  on  the  grid. 

In  the  case  of  a  square  grid  [1, . . . ,  N]d ,  K  =  N/2.  The  result  may  be  dismaying 

at  first  glance  because  it  appears  that  the  approximation  error  dh  —  bh  must 

1 

be  decreased  faster  than  0(1/N2)  just  to  remain  stable.  However,  this  type  of 
bound  is  natural,  since  the  operators’  norms  themselves  are 

iin*«  Map  =  o(k2). 

The  large  operator  norm  is  normalized  by  the  stability  condition  dt  =  0(dx2),  so 
one  can  think  of  these  bounds  in  the  update  sense  as: 

|| (F-dtQ)  -  (L-dtP) ||op  <  \\ah  —  ah\\l . 

Proof.  The  result  is  a  simple  consequence  of  Young’s  inequality:  ||/*g||2  < 
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Il/lli  h\\2-  We  have 


\\Q-p\\op=  sup  ||(Q--P)«||2 
ink 

=  ||fc(ah  -  a'h)  *  (ku) ||2 


<  \\k(ah  -  a^WkuWz 

<  K2  \\ah  -  a!^  . 


□ 

The  proof  clearly  generalizes  to  hyperbolic  operators  of  the  form  Qu  =  a*(ik  u) 
as  well. 

For  an  example,  recall  the  forward  Euler  discretization  of  a  parabolic  PDE: 

un+l  =  (/  -  dtLh)un, 

over  a  time  interval  [0,T],  If  \\dh  —  a,h\\1  =  5,  approximating  Q  by  P  incurs  an 
additional  local  truncation  error  of  magnitude  8K2dt  at  each  time  step. 

As  the  grid  is  refined,  the  CFL  condition  requires  that  K2dt  stay  approxi¬ 
mately  constant,  so  that  the  approximation  error  per  step  remains  approximately 
constant. 

2.6  The  Modified  Equation  Prospective 

Using  the  variational  principle  for  the  explicit  scheme  applied  to  the  parabolic 
equation  yields  the  following  first  order  optimality  condition: 

0  G  un+1  —  (1  —  dtLh)un  —  dtfh  +  fid  ||un+1||1 , 

which  is  equivalent  to 

/\n+l  —  n.n  ^  *  II 

0  e  — F - +  L^n  ~  h  +  ^<9  \\un+1  . 

at  dt  "  111 
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Taking  fi  =  Sdt  and  formally  sending  dt  and  h  to  zero  leads  to 

ut  +  Lu  —  /  G  —5d  Hfi ||  1 ,  (2-11) 

or 

ut  +  Lu  =  /  —  5p{u)  (2-12) 

where  p{u)  denotes  the  particular  element  of  the  subdifferential  so  that  the  differ¬ 
ential  inclusion  (2.11)  is  an  equality.  The  sparse  scheme  applied  to  hyperbolic  and 
elliptic  problems  yields  analogous  modified  equations.  We  consider  this  to  be  the 
modified  equation  in  the  sense  that  the  numerical  scheme  is  directly  solving  this 
problem.  The  subgradient  contribution  is  a  vanishing  ‘compression’  term  which 
may  be  interpreted  as  a  force  which  pushes  the  solution  u  toward  the  nearest  (in 
the  L 1  proximal  sense)  union  of  low  dimensional  subspaces  spanned  by  the  Fourier 
basis. 

Well-posedness  for  the  modified  equation  is  guaranteed  via  the  theory  of  differ¬ 
ential  inclusions  on  Banach  spaces  ( e.g .  [9,  26]).  The  theorem  below  summarizes 
these  results  in  the  current  context. 

Theorem  2.6.1  (Well-posedness).  Let  u(t)  satisfy  the  differential  inclusion 

d tu{t)  G  -A(u(t))  -  5d\\u(t)\\Li 

with  -u(O)  in  the  domain  of  the  monotone  (single-valued)  operator  A.  Then  for  all 
5  >  0  ,  there  exists  a  unique  solution  u(t)  defined  for  all  t  >  0  which  is  the  solution 
to 


dtu(t)  =  — A(u(t ))  -  Sp(u(t)) 

for  some  p  G  <9||-u(f)||Li . 

Lastly,  we  mention  that  if  we  want  to  directly  compare  the  error  between  the 
solutions  of  the  original  and  modified  equations,  the  error  grows  linearly  in  time 
(at  worst). 
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Theorem  2.6.2.  Let  u  be  the  solution  to 


and  let  us  solve 

Then 


ut  +  Lu  =  f 


(■ us)t  +  Lus  -  f  G  5d  || m,5  || ! 

II u(t,  •)  -  us{t,  •) 1 1 2  <  2 St. 


The  proof  is  direct  and  can  be  found  in  the  Appendix.  Similar  results  are 
easily  proved  for  the  elliptic  and  hyperbolic  cases  using  only  that  ||p(tt)||oc  <  1 
and  standard  energy  estimates;  this  approach  also  provides  a  simple  alternate 
proof. 

2.7  Denoising  Perspective 

Soft  thresholding  also  appears  in  early  methods  for  signal  denoising  using  wavelets 
[30].  We  refer  the  reader  to  that  work  for  full  details,  and  list  here  only  the 
analogues  of  its  major  results  in  the  current  context. 

Consider  the  following  denoising  problem:  we  wish  to  recover  a  signal  /  G 
Mn  from  noisy  observations  d  =  f  +  w,  ||ie||i  <  /i,  by  soft-thresholding  DFT 
coefficients  by  /j.  This  approach  enjoys  the  following  properties: 

•  (Smoothing)  The  recovered  signal  flt  satisfies  1 1 /M 1 1 Hk  <  ||/||Hfc  for  any 
Sobolev  norm  ||-||Hfe.  In  particular,  |//i(/c)|  <  \f(k)\  for  all  frequencies  k. 

•  (Near  optimality)  ./),  is  near-minimax: 

sup  sup  Wfn-fWp  <4 inf  sup  sup  1 1 /(d)  -  / \\% 
}\f\\Hk<Ci  ||™||i<m  n  /  \\f\\Hk<Ci  |H|i<m 

where  /(d)  is  any  other  estimator  of  /. 
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The  smoothing  property  guarantees  that  the  recovered  signal  is  ‘noise-free’;  the 
near  optimality  property  guarantees  that  for  worst-case  signals  of  bounded  Sobolev 
norm  and  noise  of  bounded  i1  norm,  the  result  recovered  by  soft  thresholding  is 
nearly  the  ‘optimal’  (see  [30]). 

Next,  consider  the  solution  u€  to  the  standard  parabolic  multiscale  problem 

=  0  on  [0,  27t]  periodic,  ue(x,  0)  =  ue0(x). 

The  theory  of  asymptotic  homogenization  {e.g.  [70])  can  be  used  to  show  that  at 
time  point  tn,  the  exact  solution  ue  satisfies 

u€(x,tn)  =  u0(x,tn)  +  eui(x,x/e,tn)  +  e2R(x,tn ) 

with  \R(x,  t)  <  C.  This  expansion  is  valid  as  long  as  we  assume  that  the  equation 
is  taken  on  a  periodic  domain  and  a{x )  is  as  smooth  as  we  like.  For  a  numerical 
solution,  the  asymptotic  expansion  can  be  easily  modified  to  include  truncation 
error  rn+1  as  follows:  if  we  let  vn+l  denote  the  numerical  solution  at  time  tn+1, 
then 

vn+1  =  uo(x,  tn+l )  +  eui(x,  x/e,  tn+l)  +  e2R(x,  tn+1 )  -  rn+1. 

This  form  allows  us  to  draw  a  connection  between  the  denoising  and  homoge¬ 
nization  problems:  for  an  appropriate  threshold  choice  //,  the  compressive  spectral 
method  denoises  vn+1  as 

vn+1  =  u0(x,  fn+1)  +  eUl(x,  x/e,  tn+l)  +  e2R(x,  tn+1 )  -  rn+1 

N - V - '  s - V - ^ 

signal  noise 

and  attempts  to  recover  the  first  terms  in  the  asymptotic  expansion.  These  inter¬ 
pretation  is  valid  between  any  two  time  steps,  but  may  not  hold  globally. 

2.8  Efficient  Implementation 

In  this  section,  we  describe  important  details  pertaining  to  the  numerical  method 
and  algorithm  considerations.  Using  a  concrete  example,  we  show  that  a  favorable 


d ue  d  (  .  due 

~  di{a{x’x/e)te 
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complexity  can  be  achieved. 


2.8.1  The  Proximal- Galerkin  Algorithm 

The  implicit  scheme  described  above  requires  fast  minimization  of  the  energy 
(2.8),  and  differs  from  many  cases  where  L 1  regularization  is  added  because  the 
problem,  e.g.  compressed  sensing  [17],  TV  minimization  [75],  or  basis  pursuit  [21], 
is  ill-posed  without  it.  For  the  multiscale  PDE  problem,  this  is  not  the  case  since 
an  appropriately  discretized  version  of  (2.8)  will  be  well-posed  and  can  be  solved 
by  inverting  a  linear  system 

Qu  =  f 

where  Q  is  a  positive-definite  (and  even  sparse,  in  physical  rather  than  Fourier 
space)  matrix.  If  the  elliptic  operator  is  discretized  appropriately,  fast  and  exten¬ 
sively  studied  preconditioned  conjugate  gradient  solvers  are  available.  So,  to  be 
competitive,  the  compressive  implicit  scheme  must  leverage  sparsity  of  the  solu¬ 
tion  u  to  perform  the  (approximate)  linear  inversion  Qu  =  /  quickly.  For  this 
purpose,  we  propose  the  hybrid  proximal  gradient  descent  and  Galerkin  approxi¬ 
mation  algorithm  described  below,  which  is  related  to  the  procedure  described  in 
[28] . 

First,  let  D  be  the  diagonal  part  of  Q.  Since  0  is  the  matrix  corresponding  to 
a  Fourier-space  discretized  elliptic  operator,  D  is  the  matrix  corresponding  to  a 
multiple  the  Fourier-space  discretized  Laplacian.  We  take  n  ~  10,  n  >  0,  u  >  0, 
and  initialize  the  solution  to  be  zero  (he.  u  —  0). 

The  Proximal-Galerkin  Algorithm 
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for  j  =  l:  n  do 

u  =  shrink (u  +  oj/1”1(/  —  Qu),  fi) ; 


end  for 

set  /  =  supp(w); 
set  u  =  arg  min  w  : 
Return  u. 


supp(w)  C  I]- 


Qw  -  f 


2 

2 


The  algorithm  begins  with  a  few  iterations  of  the  proximal  gradient  method 
applied  to  the  energy 

7-1/  \  II  'll  I  1  lT l  IT  £l 

b{w)  =  n\\w  \\1  +  -w  Qw  —w  j 

where 

Q'  =  D~1/2QD~1/2, 

f  =  D-X'2h 

w'  =  D1/2w. 

This  is  a  simple  Jacobi  preconditioning  of  the  analogous  energy  with  Q,  w ,  and 
/.  Rather  than  iterating  proximal  gradient  to  convergence,  which  would  be  too 
slow,  the  algorithm  stops  after  just  a  few  iterations  with  rough  approximation. 
The  support  of  that  solution  is  used  to  identify  the  Fourier  modes  with  largest 
magnitude  coefficients,  and  then  a  Galerkin  approximation  is  computed  over  those 
modes.  Due  to  sparsity  in  the  Fourier  domain,  the  linear  solve  associated  with 
the  Galerkin  part  is  small  and  inexpensive-  computational  complexity  depends 
on  the  grid  mesh  size  only  through  the  sparsity  of  the  solution. 

2.8.2  Algorithm  Complexity 

The  pseudospectral  approach  of  computing  the  convolution 

ka*  (ku) 
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uses  an  FFT,  and  for  an  iV-gridpoint  problem  this  reduces  the  computational 
complexity  per  iteration  from  0(N2)  to  0(N  log  N).  We  now  consider  the  com¬ 
putational  complexity  of  the  sparse  spectral  method,  which  must  be  comparable 
to  O(iVlogiV)  to  be  practical. 

Suppose  that  the  sparsely  approximated  operator  is  defined  Pu  —  ka!  *  (ku), 
where  the  sparsity  (number  of  nonzeros)  of  a '  is  m,  and  that  the  sparsity  of  u  is 
r.  By  treating  the  a'  *u  sparse  convolution  as  a  summation  of  sparse  vectors,  it 
can  be  accomplished  with  complexity 

0(mr  min(logr,  logm)),  (2-13) 

free  of  any  dependence  on  the  full  problem  size  N,  by  storing  the  sparse  vectors 
a!  and  u  as  sorted  linked  lists  and  computing  the  sum  as  a  merge  operation,  with 
a  priority  queue.  For  the  modest  one-time  cost  of  initializing  a  length  N  array, 
the  complexity  can  be  decreased  to 

0(mr)  (2-14) 

by  leaving  the  sparse  vectors  unsorted.  We  iterate  over  the  mr  nonzero  coefficients 
which  must  be  added,  and  use  an  auxiliary  array  keep  track  of  the  partial  result. 
When  a  new  coefficient  of  the  partial  result  becomes  nonzero,  it  is  placed  in  a 
growing  list  of  indices.  After  we  have  visited  each  of  the  mr  coefficients  to  be 
added,  we  iterate  over  the  list  of  nonzero  indices,  perform  the  shrink  operation  on 
the  corresponding  auxiliary  array  entry  holding  the  partial  result,  and  copy  the 
outcome  into  a  list  which  holds  the  final  result.  Along  the  way,  we  ‘zero  out’  the 
entry  of  the  partial  result  array,  never  incurring  another  O(N)  cost. 

Finally,  if  the  problem  is  elliptic  or  requires  implicit  time  steps  and  the  Proximal- 
Galerkin  algorithm  is  used,  the  complexity  includes  a  term 

0(r 3), 

the  cost  of  the  Galerkin  linear  solve  over  the  support  found  with  proximal  gradient. 
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Both  (2.13)  and  (2.14)  are  preferable  to  the  0(N  log  N)  cost  of  the  pseudospec- 
tral  method  for  very  sparse  problems  and  in  the  homogenization  limit  discussed 
next  in  Section  2.8.3.  For  the  numerical  examples  considered  in  this  limit,  m 
and  r  stay  approximately  constant,  leading  to  computation  time  which  does  not 
increase  as  the  grid  is  refined. 

One  key  to  the  effective  application  of  the  sparse  spectral  method  is  proper 
discretization.  For  a  typical  homogenization  problem,  we  are  interested  in  the 
solution  of  an  equation  such  as 


— div(a(x/e)  Vu)  =  / 

for  e  close  to  zero,  and  we  might  choose  the  inhomogeneity  coefficient 

1 

a(x )  =  1  +  -  sin  nx. 

This  choice  is  ideally  sparse  in  the  Fourier  domain,  with  only  three  nonzero  entires 
regardless  of  N ,  using  the  standard  uniform  grid.  If  e  =  1/1000,  then  a  still  has 
only  three  non  zeros.  However,  choosing  e  =  707^  results  in  a  being  completely 
dense.  These  two  choices  of  e  differ  by  less  than  10~6,  and  the  first  leverages 
extreme  sparsity  in  the  problem  while  the  second  does  not.  This  example  shows 
that  it  is  prudent  to  assume  a  certain  relationship  between  the  grid  spacing  and 
e,  considered  next. 

2.8.3  Homogenization  Limit 

For  homogenization  problems  in  particular,  where  one  is  interested  in  the  limit 
e  — y  0,  we  can  keep  Ne  fixed  as  the  grid  is  refined.  Empirically,  we  have  observed 
that  this  keeps  the  sparsity  of  the  operator  and  of  the  solution  approximately 
constant.  For  a  simple  case  of  this  Ne  =  c  (c  a  constant)  limit,  the  following 
theorem  guarantees  the  sparsity  of  the  operator  remains  fixed  along  a  subsequence. 


75 


Theorem  2.8.1.  Let  Le  be  the  elliptic  operator  defined 


Lev  =  —  div(a(x/e)Vu), 


and  let 


Q€,nu  =  kafi  *  (ku) 

be  its  Fourier  discretization  on  an  N— point  discretization  of  [0,  2tt).  Then  Q6,n 
and  Qe/2,2N  «re  equally  sparse:  that  is, 


#{k  :  \a,2N(k)\  >  A}  =  #{/c  :  \aN(k)\  >  A/2}  (2.15) 


for  all  A  >  0. 


See  the  appendix  for  a  proof.  Note  that  the  theorem  assumes  the  standard 
definition  of  the  DFT  on  N  grid  points, 

N-i  ,  . 

J7v[a(a;)](/c)  = 

j= o  ^ 

which  is  not  unitary.  This  accounts  for  the  appearance  of  A/2  rather  than  A  on 
the  right  hand  side  of  (2.15).  This  factor  cancels  out  in  the  end  because  with  this 
definition  of  the  DFT,  the  il  norm  in  Theorems  2.5.5  and  2.5.6  should  be  scaled 
by  1/Ah 

The  complexities  (2.13)  and  (2.14)  become  very  favorable  in  the  Ne  =  c  limit, 
where  m  and  r  remain  nearly  constant  or  grow  approximately  logarithmically 
with  N  as  the  grid  is  refined.  In  each  case  we  observed,  the  overall  algorithm 
complexity  is  linear  or  sub-linear  in  N. 


^  g-27 njk/N 


2.9  Numerical  Examples 

In  [76],  the  authors  demonstrated  the  effective  application  of  the  compressive 
spectral  method  to  a  variety  of  problems.  Here,  we  expand  on  those  results  and 
give  examples  of  the  additions  to  the  method  proposed  in  this  chapter:  the  implicit 
scheme  and  sparse  operator  approximation. 
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2.9.1  Transport  Equation,  ID 


The  PDE  considered  is  the  traveling  wave  equation: 


ut  +  a(x)ux  =  0, 
x  G  [0,  27r]  periodic, 
u(x,  0)  =  sin(x) 


with  oscillatory  coefficient 


.  .  1  (  0.6 +  0.2  cos  a: 

a[x)  =  -  exp 


1  +  0.7  sin  128a: 
The  update  is  given  by  leap  frog  time  discretization: 


u 


=  un  1  —  2 dtta  *  ( ikun ). 


We  choose  the  above  form  for  a  throughout  this  section,  because  it  is  less  sparse 
than  simple  trigonometric  functions. 

The  grid  sizes  considered  are  N  =  210, . . . ,  214  and  the  values  of  other  param¬ 
eters  are  dtt  =  6.25 x  10-6,  ||a  —  a'l^  =  10-2,  n  =  1.2 x  10-5,  and  the  simulation  is 
run  to  a  final  time  t  =  0.5. 

Figure  2.3  shows  the  full  spectral  and  compressive  spectral  (sparse  opera¬ 
tor/sparse  solution)  solutions  on  coarse  and  fine  scales.  The  compressive  scheme 
correctly  captures  the  largest  Fourier  coefficients  of  the  solution,  discarding  all  but 
3.7%,  and  the  operator  approximation  discards  all  but  2.6%.  The  “true"  solution 
was  computed  on  a  fine  grid  with  finite  difference  methods. 

Figure  2.4  shows  the  L2  error  and  sparsity  of  the  compressive  spectral  approx¬ 
imations  as  the  grid  is  refined  with  dt  held  constant.  Error  is  computed  as  the 
L2  distance  to  the  full  spectral  solution.  The  error  of  the  sparse  operator/sparse 
solution  scheme  is  dominated  by  the  sparse  approximation  of  the  solution;  spu¬ 
rious  modes  in  the  leap  frog  scheme  make  a  sparse  approximation  of  it  difficult. 
Over  the  range  of  grids  considered,  sparsity  of  the  operator  eventually  becomes 
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constant  while  sparsity  of  the  solution  grows  about  linearly.  The  complexity  of 
the  compressive  spectral  method  is  thus  linear  in  N  over  the  grid  sizes  considered. 

Figure  2.5  considers  the  same  problem  but  with  a  resonant  forcing  term 

f(x)  =  esin{x/128)2 

with  N  =  2048  and  all  other  parameters  the  same  as  the  non-forced  problem.  The 
solution  has  11.3%  nonzero  Fourier  coefficients,  with  ||ufun  —  uSparse||2  =  2.5xl0-3. 
The  resonant  forcing  causes  sharp  and  irregular  oscillations  at  the  fine  scale,  which 
make  the  problem  less  sparse,  but  the  compressive  scheme  still  captures  the  correct 
behavior. 


2.9.2  Elliptic  Problem,  ID 


The  PDE  considered  is  the  elliptic  problem: 

—  (a(x)ux)x  =  sin  2a;, 
x  €  [0,  27 r]  periodic, 
/  udx  =  0 


with 


a(x)  =  exp 


0.6  +  0.2  cos  a; 


1  +  0.7  sinx/e/ 

such  that  Ne  =  8,  and  usual  spectral  operator  discretization: 


Lhu  =  kah*  (ku)  =  /. 


This  time  we  consider  the  homogenization  limit,  keeping  Ne  =  8  with  e  = 
6+  I+b  •••  >  hid’  and  se^  ||d  —  a'||i  =  lx  10-4.  Parameter  values  for  the  Proximal- 
Galerkin  algorithm  are  n  —  10,  n  —  5x  10-8,  and  w  =  5x  10-3. 

Figure  2.6  shows  the  full  spectral  and  compressive  spectral  solutions  on  coarse 
and  fine  scales.  Both  the  sparse  solution  and  operator  approximation  keep  8.5% 


78 


k 

Figure  2.3:  Left:  True  (blue)  and  sparse  operator/sparse  solution  (green)  solu¬ 
tions  in  physical  space.  The  two  curves  lie  almost  on  top  of  each  other.  Right: 
Zoomed  in  true  (blue)  and  sparse  (green  ‘x’)  solutions.  Bottom:  True  (blue)  and 
sparse  (red  ‘o’)  solutions  in  Fourier  space.  N  =  4096,  operator  nonzeros  =  107, 
solution  nonzeros  =  153. 
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Figure  2.4:  Left:  Sparse  operator/full  solution  (blue),  full  operator/sparse  so¬ 
lution  (green,  dashed),  and  sparse  operator/sparse  solution  (red  x)  L2  distance 
to  the  full  spectral  solution  as  the  grid  is  refined.  The  y  axis  has  a  log10  scale. 
Right:  Number  of  nonzero  Fourier  coefficients  of  the  operator  (blue)  and  solution 
(green,  dashed)  as  the  grid  is  refined.  The  y  axis  has  a  log2  scale. 

of  the  coefficients.  Note  that  the  full  result  and  sparse  operator /sparse  solution 
result  lie  almost  on  top  of  each  other,  even  at  the  resolution  of  the  fine  scale. 

Figure  2.7  shows  error  ( L 2  distance  to  the  full  spectral  solution)  and  sparsity 
under  refinement.  “Sparse  operator"  refers  to  the  solution  obtained  with  the 
sparsely  approximated  operator,  using  either  a  high  accuracy  conjugate  gradient 
solve  or  the  Proximal- Galer kin  algorithm.  “Sparse  solution"  refers  to  the  use  of 
the  Proximal-Galerkin  algorithm,  with  either  the  full  or  sparse  operator. 

Approximation  error  does  not  increase  while  both  solution  and  operator  spar¬ 
sity  remain  approximately  constant,  leading  to  computation  time  approximately 
independent  of  N.  With  N  =  213,  the  sparse  approximation  maintains  six  digits  of 
accuracy  with  only  1.1%  of  the  coefficients  of  both  the  operator  and  the  solution. 

Figure  2.8  illustrates  that  for  a  fixed  number  of  nonzero  coefficients,  the  sparse 
operator  approximation  incurs  smaller  error  than  the  solution  approximation. 
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Figure  2.5:  Left:  True  (blue)  and  sparse  operator/sparse  solution  (green)  solu¬ 
tions  with  resonant  forcing  term  in  physical  space.  Right:  Zoomed  in  true  (blue) 
and  sparse  (green  ‘x’)  solutions.  Bottom:  True  (blue)  and  sparse  (red  ‘o’)  solu¬ 
tions  in  Fourier  space.  N  =  2048,  operator  nonzeros  =  86,  solution  nonzeros  = 
231. 
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Figure  2.6:  Left:  True  (blue)  and  sparse  operator/sparse  solution  (green)  solu¬ 
tions  in  physical  space.  The  small  rectangle  shows  the  axis  limits  of  the  zoomed 
in  plot  to  the  right.  Right:  Zoomed  in  true  (blue)  and  sparse  (green  ‘x’)  so¬ 
lutions.  Bottom:  True  (blue)  and  sparse  (red  ‘o’)  solutions  in  Fourier  space. 
N  =  1024,  operator  nonzeros  =  86,  solution  nonzeros  =  87. 
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Figure  2.7:  Left:  Sparse  operator/full  solution  (blue),  full  operator/sparse  solu¬ 
tion  (green,  dashed),  and  sparse  operator/sparse  solution  (red  x)  error  under  the 
homogenization  limit.  The  y  axis  has  a  log10  scale.  Right:  Number  of  nonzero 
Fourier  coefficients  of  the  operator  (blue)  and  solution  (green,  dashed)  as  the  grid 
is  refined. 


log2(gridpoints) 


Figure  2.8:  Pareto  curves  showing  the  tradeoff  between  approximation  error  and 
sparsity  of  the  operator  (blue)  and  solution  (green,  dashed). 
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2.9.3  Parabolic  Problem,  ID 

The  PDE  we  consider  here  is  the  parabolic  equation: 

ih  ~  (d(x) ii =  0, 
x  G  [0,  27r]  periodic, 
u(x,  0)  =  1  +  cos(a:  —  n) 

with  /0.6  +  0.2  cosx\ 

a  (x)  =  exp  - — 

\1  +  0.7  smx/e  J 

We  again  consider  the  Ne  =  8  limit,  e  =  and  set  ||a  —  a'l^  = 

1  x  10-2  and  dt  —  lx  10~2  for  all  N.  Parameter  values  for  the  Proximal-Galerkin 
algorithm  are  n  =  10,  //  ranges  from  5xl0-6  to  6.4  xlO”6,  and  u  =  lxlO-2. 

Figure  2.9  compares  the  solutions  on  coarse  and  fine  scales.  The  sparse  solu¬ 
tion  retains  3.2%  of  the  coefficients  and  the  operator  is  also  approximated  with 
3.2%.  Figure  2.10  shows  error  and  sparsity  under  refinement.  Approximation 
error  decreases  while  sparsity  of  both  operator  and  solution  stay  constant.  The 
overall  complexity  is  thus  constant  in  N  over  the  range  of  grid  sizes  considered. 
For  this  problem,  sparse  approximation  of  the  operator  incurs  most  of  the  error. 


2.9.4  Elliptic  Problem,  2D 


We  consider  the  elliptic  problem 


-div(a(x)Vu)  =  10  sin  x  sin?/, 
x,  y  G  [0,  27t]  periodic, 

[  u  dxdy  =  0 


with 


,0.6  +  0.2  cos  x  0.6  +  0.2  cos  y 
a(x’v)  =  6XP  I  1  +  0.7 sin  r/e  +  l  +  O.Tsmy/e 
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Figure  2.9:  Left:  True  (blue)  and  sparse  operator/sparse  solution  (green)  solu¬ 
tions  in  physical  space.  Right:  Zoomed  in  true  (blue)  and  sparse  (green  ‘x’) 
solutions.  Bottom:  True  (blue)  and  sparse  (red  ‘o’)  solutions  in  Fourier  space. 
N  =  2048,  operator  nonzeros  =  64,  solution  nonzeros  =  65. 
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Figure  2.10:  Left:  Approximation  error  of  the  sparse  operator/full  solution  (blue), 
full  operator/sparse  solution  (green,  dashed),  and  sparse  operator/sparse  solution 
(red  x)  error  under  the  homogenization  limit.  The  y  axis  has  a  log10  scale.  Right: 
Number  of  nonzero  Fourier  coefficients  of  the  operator  (blue)  and  solution  (green, 
dashed)  are  constant  as  the  grid  is  refined. 

on  an  N  x  N  grid  such  that  Ne  =  8,  with  e  =  ... ,  ^  and  ||a  —  a'l^  =  1. 

Parameter  values  for  the  Proximal-Galerkin  algorithm  are  n  =  20,  /j  between 
4 x  1(T4  and  32 x  1(T4,  and  w  =  2x  10"2. 

Because  the  full  and  spectral  solutions  are  very  close  to  each  other  in  phys¬ 
ical  space  and  an  overlaid  comparison  of  surfaces  is  difficult,  Figure  2.11  shows 
the  solutions  on  a  log  scale  in  Fourier  space.  Of  the  220  coefficients  in  the  full 
solution,  the  sparse  solution  and  operator  retain  just  0.2%  while  maintaining  four 
digits  of  accuracy.  Figure  2.12  shows  that  approximation  error  decreases  slightly 
with  constant  sparsity  and  computation  time.  For  some  grid  sizes,  the  sparse 
operator/sparse  solution  scheme  actually  attains  a  lower  error  than  the  sparse  op¬ 
erator/full  solution  scheme,  evidence  of  the  denoising  effect  discussed  in  Section 
2.7. 

To  compare  the  Fourier  coefficients  of  the  full  and  sparse  solutions  more  ac- 
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Figure  2.11:  Full  (left)  and  sparse  (right)  solutions  on  a  log  scale  in  Fourier  space. 
Note  that  the  great  majority  of  coefficients  in  the  sparse  solution  are  exactly  zero. 
N  =  1024,  e  =  p^g,  operator  nonzeros  =  1972,  solution  nonzeros  =  1874. 

curately,  the  left  panel  of  Figure  2.13  shows  the  magnitude  of  the  4500  largest 
Fourier  coefficients  of  the  true  solution  sorted  in  descending  order.  The  magni¬ 
tude  of  the  corresponding  sparse  solution  Fourier  coefficients  is  also  shown,  with 
an  upward  bias  to  account  for  all  the  wave  numbers  not  present.  The  right  panel 
shows  the  fraction  of  full  solution  wave  numbers  which  are  captured  by  the  sparse 
scheme.  The  compressive  scheme  correctly  identifies  all  500  of  the  largest  modes 
in  the  full  solution,  and  about  68%  of  the  full  solution’s  largest  1800  modes. 


2.10  Conclusion 

In  this  chapter,  we  have  proposed  a  sparse  operator  approximation  and  an  effi¬ 
cient  method  for  extending  the  work  of  [76]  to  implicit  solvers  (Section  2.4).  We 
have  proven  the  convergence  of  the  original  compressive  spectral  scheme  [76]  and 
the  new  variants,  including  a  modified  equation  which  shows  the  effect  of  soft 
thresholding  is  equivalent  to  including  an  L 1  subgradient  term  in  the  PDE.  Also, 
we  connect  the  homogenization  problem  with  that  of  signal  denoising  via  wavelet 
thresholding.  For  PDE  with  sparse  initial  data  or  forcing  terms,  the  new  methods 
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Figure  2.12:  Left:  Approximation  error  of  the  sparse  operator/full  solution  (blue), 
full  operator/sparse  solution  (green,  dashed),  and  sparse  operator/sparse  solution 
(red  x)  error  under  the  homogenization  limit.  The  y  axis  has  a  log10  scale.  Right: 
Number  of  nonzero  Fourier  coefficients  of  the  operator  (blue)  and  solution  (green, 
dashed)  are  constant  as  the  grid  is  refined. 


Figure  2.13:  Left:  energy  spectrum  decay  of  the  full  and  sparse  solutions.  The 
plot  shows  just  the  largest  4500  coefficients  of  the  full  solution,  the  support  of 
which  contains  all  coefficients  of  the  sparse  solution.  Right:  fraction  of  sparse 
modes  appearing  among  the  largest  n  true  modes,  as  a  function  of  n. 
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are  asymptotically  preferable  to  the  pseudospectral  approach.  The  methodol¬ 
ogy  presented  here  could  be  translated  to  other  psuedospectral  methods  which 
employ  alternative  bases.  Computationally,  this  amounts  to  replacing  the  Fast 
Fourier  transforms  in  the  psuedo-codes  above  with  the  appropriate  transforma¬ 
tion.  This  could  be  useful  in  cases  where  the  solutions  are  sparse  against  another 
known  basis. 

2.11  Appendix 

Before  giving  the  proofs  of  the  theorems  from  Section  2.5.1,  we  recall  the  definition 
of  Bregman  Distance  (also  known  as  Bregman  Divergence). 

Definition  2.11.1.  Let  J  be  a  convex  function  and  u,v  be  points  in  the  domain 
of  J.  Also  let  p  be  an  element  of  the  subdifferential  of  J,  i.e.  p  G  d J(v).  We 
define  the  Bregman  Distance  between  u  and  v  as 

Dj(u,  v )  =  J(u )  —  J(v )  —  (p,u  —  v ). 

In  general,  the  Bregman  Distance  is  not  symmetric  and  does  not  obey  the 
triangle  inequality,  so  it  is  not  a  distance  in  the  typical  sense. 

In  what  follows,  we  will  also  use  basic  facts  regarding  monotone  operators. 

Definition  2.11.2.  Let  A  be  a  multi-valued  map  from  V  into  itself.  We  call  A 
monotone  if  and  only  if  for  any  u,v  G  Dom(A)  and  any  values  Au  and  Av  might 
take  on, 

(u  —  v,  Au  —  Av)  >  0. 

If  A  =  d F  is  the  subdifferential  of  a  convex  function,  then  it  is  monotone. 


We  can  now  give  the  proof  of  theorem  2.5.1,  in  which  we  omit  hats  for  nota- 
tional  clarity. 
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Proof.  Consider  the  iterations  for  arbitrary  points  un  and  vn: 

7/n+ 1  _  n,n 

M«n+')  +  — =  -L„un  +  f 

v n+1  _  n,n 

ra>K+1)  +  — s +  /• 

By  taking  the  difference  between  these  two  equations  we  arrive  at 

p(p(un+1)  -p(vn+1))  +  ^(«n+1  -  ^n+1)  -  Jt(un  -  vn)  =  —Lh{un  -  vn) 

and  taking  the  inner  product  of  this  equation  with  un+1  —  vn+1  yields 

H  ( p(un+1 )  -  p(vn+1),  un+1  -  vn+1)  +  -  (un+1  -  vn+1,  un+1  -  vn+1)  - 
'  'at  '  ' 

J  < un  -  vn ,  un+1  -  vn+1)  =  l\—Lh{un  -  vn),un+1  -  vn+1"j  . 

Rearranging  terms  and  taking  upper  bounds  we  get  the  following: 

Hdt  (p(un+1)  ~p(vn+1),un+1  -vn+1)  +  \\un+1  -  vn+1\\2 

=  < un  -  vn,un+1  -  vn+1)  +  (~dtLh(un  -  vn),un+1  -  vn+ ^ 

=  ((I-  dtLh)(un  -  vn)lUn+l  -  vn+1^ 

<  1 1  (/  -  dt Lh)(un  -  vn)  1 1  j  | un+1  -  vn+1 1 1 

<  1 1 (/  -  dtLh)\\op\\un  -un|| ||un+1  -  vn+1 1 1 . 

Note  that  pdt  (p(un+1 )  —  p(vn+1) ,  un+1  —  vn+l)  is  non-negative  by  monotonicity  of 
the  subgradient  of  a  convex  function.  We  show  this  here  by  using  the  nonnegativity 
of  Bregman  distance: 

0  <  DpF(un+1,vn+1)  +  DpF(vn+1,un+1) 

=  F(un+l )  -  F(vn+1)  -  (p(vn+1),un+1  -  vn+1 >  +  F(vn+1)  -  F(un+1)  -  (p(un+1) ,  vn+1  -  un+1 ) 
=  < p{un+l )  -  p(vn+1),un+1  -  vn+1)  . 

Combining  the  positivity  of  the  subgradient  terms  with  equation  above  provides 
us  with  the  following  bound  (assuming  || (/  —  dtLh) \\op  <  1): 

\\un+1  -  vn+1\\  <  \\{I  -  dtLh)\\op\\un  -  vn\\  <  \\un-vn\\ 

as  desired.  □ 
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Proof  of  theorem  2.5.2. 


Proof.  Considering  the  optimality  condition  for  the  energy  (2.8)  defining  the  im¬ 
plicit  scheme,  we  see  that  the  iterations  for  un  and  vn  can  be  written 

n  .ra  +  1  _  n.n 

pp{un+l)  + - - - =  -LhUn+1  +  f 

7,n+ 1  _  n.n 

w(vn+1)  +  — —  =  -Lhvn+1  +  f . 

(If  the  operator  Lh  being  considered  in  (2.8)  is  not  positive  semidehnite,  then  use 
(2.7)  instead.)  By  taking  the  difference  between  these  two  equations  we  arrive  at 

li(p(un+1)  -  p(vn+1))  +  Jt(un+l  -  vn+1)  -  Jt(un  -  vn)  =  —Lh{un+1  -  vn+l). 

Next,  taking  the  inner  product  of  this  equation  with  un+ 1  —  vn+1  yields 

II  ( p{un+1 )  -  p(vn+1),un+1  -  vn+1)  +  (un+1  -  vn+1,  un+1  -  vn+1)  - 

LLL 

- 1 

—  < un  -  vn ,  un+1  -  vn+1)  =  (-Lh(un+1  -  vn+1),  un+1  -  vn+l)j 

Re-arranging  terms  and  taking  upper  bounds  we  get  the  following: 

(idt  ( p{un+l )  -p(vn+1),un+1  -  vn+1 )  +  || un+1  -  un+1||2 
=  (un  -  vn,  un+1  -  vn+1)  +  (-dtLh(un+1  -  vn+1),  un+1  -  vn+lXj  . 

As  in  the  explicit  timestep  case,  (p(un+1)  —  p(vn+1),  un+l  —  vn+l)  >  0  and  so 

|| un+1  -  vn+l\\2  <  (un  -  vn,un+1  -  vn+1)  +  (-dtLh{un+l  -  vn+1),un+1  -  vn+1^  . 

If  Lh  is  positive  semidehnite  then  we  have 

||un+1  -un+l||2  <  (un  —  vn,  un+1  —  vn+1)  <  ||un  -un||||un+1  -un+1|| 

and  by  canceling  out  terms  we  get  the  contractive  inequality 

||«n+1  -Un+1||  <  \\un-vnW 

as  desired.  □ 
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Proof  of  Theorem  2.5.3: 


Proof.  We  assume  that  S  is  stable  in  the  following  sense: 

||wn+1||  <  ||hn|| 

for  some  lp  norm;  common  choices  would  be  the  l 2  or  Z°°  norms.  Because  the 
shrink  operator  decreases  the  magnitude  of  each  component  of  a  vector,  it  will 
(strictly,  because  //  >  0)  decrease  whatever  norm  is  chosen  (in  fact,  the  shrink 
operator  is  a  contraction  in  all  lp  norms),  ft  follows  easily  that 

',rc+1"  <  \\Q(un  . ..  ,ur. 


u 


<  \\un\ 

fl  II  —  II  ^  •  •  •  5  W/X  /II  —  lljul 


"n—k\ 

lH  ) 


so  that  the  stability  of  S  implies  the  stability  of  S M.  In  fact  StJ  is  more  stable  than 

S. 

The  key  observation  for  showing  consistency  of  Su  is  that  while  shrink(-,/i)  is 
nonlinear,  the  amount  of  this  nonlinearity  is  bounded.  In  particular, 


shrink(x,  /i)  =  x  +  0(/i) 

for  any  x,  with  |0(/u)j  <  fi.  Applying  this  observation  to  the  definition  of  the 
sparse  scheme  and  assuming  (for  the  purpose  of  local  truncation  analysis)  that 
both  schemes  have  the  same  starting  points  h”  =  un, . . . ,  u'^pk  =  un~k, 

^+1  =  shrink(g(^,...,^-fc),/i) 

=  Q(u;,...,un-k)  +  0^) 

=  Q(un,...,un-k)  +  0(fi) 

=  hn+1  +  0(/i). 

This  shows  that  locally,  S  and  S f  differ  only  by  a  0(/i)  quantity.  This  quantity 
may  naively  be  accounted  as  part  of  the  local  truncation  error  for  the  sparse 
scheme,  in  which  case 

r;  =  T"  +  0(/i) 
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where  rn  denotes  the  local  truncation  error  of  S  and  r,?  the  local  truncation  error 

fj, 

of  Slt. 

For  the  consistency  of  S^,  we  need  the  local  truncation  error  to  be  greater 
than  first  order;  assuming  the  consistency  of  S  and  that  //  =  0(dt1+s )  yields  this 
result.  When  /i  =  0(dtp )  for  some  p  such  that  rn  =  0(dtp)  as  well,  r”  =  0(dtp) 
and  the  order  of  convergence  of  the  scheme  is  unchanged.  □ 


Proof  of  Theorem  2.5.4: 


Proof.  First,  recall  that  the  optimality  condition  for  (2.8)  is 
dP(K+1)  +  (J  +  dtLh)un+1  -  un  +  dtfh  =  0 


".n-\- 1 1 


(2.16) 


where  p(h))+1)  €  9||w^+1||  .  For  simplicity  of  notation,  let  w  :=  (h))+1  —  -un+1). 
Assuming  (again  for  the  purpose  of  local  truncation  analysis)  that  both  schemes 
have  the  same  starting  point  u”  =  itn,  subtracting  the  ordinary  backward  Euler 
update  from  this  gives 


(/  +  dtLh)w  = 


n+ 1  \ 


which  implies 


(/  +  dtLh)w 


<  fi. 


Then,  using  the  fact  that  Lh  is  positive  definite,  we  get 


(/  +  dtLh)w 


Ml2  /n1/2 


> 


(/  +  dtLh)w 


>  1 


w 


which  gives 


\w\ 


\w\ 


L2(Q)  ~  A^1/2 


<  h- 


So,  as  with  the  explicit  scheme, 


u 


n+ 1  _  „-,n+ 1 


=  hn+i  +  O(fi)  (in  L2(Q)) 


rf  =  rn  +  O(p) 
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which  yields  consistency  if  /i  =  0(dt1+s )  with  <5  >  0,  and  implies  the  order  of 
convergence  is  the  same  as  that  of  the  ordinary  spectral  scheme  if  p,  —  0(dtp ) 
with  p  such  that  rn  =  0(dtp). 

To  prove  stability  of  the  scheme,  return  to  (2.16)  with  /  =  0  and  take  the 
inner  product  with  u™+1  to  get 

f*  HC'ili  +  + dtLh)up'  ~  (up'fu; = o 

which  leads  to 


and 


as  desired. 


u. 


<  K+1’K)  -  »  liClli  -  dt(ul+1)T Lhu\ 


<  u„ 


n+l||  -n 

2  M  2 


U 


n+1 1 


„  <  , 
2  M  2 


□ 


Proof  of  theorem  2.6.2: 


Proof.  We  have 

“ IKV)  —  us(t,  -)l ll  =  {u  -  Us,  dtu  -  dtus) 

=  (u  -  u5, -Lu  +  f  -  {-Lus  +  /  -  50||w(t)||i)) 
=  -  (u-  us,  Lu  -  +  5  (u  -  U5,50||u5||i) 

<5{u-  us,d0\\us\\i) 

<  5\\u  —  ue||2. 


It  follows  that 


—  ||u(t,  •)  ue(t,  *) 1 1 2  <  2h 


from  which  the  result  follows. 


□ 
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Proof  of  theorem  2.8.1: 


Proof.  Let  T ^[a{x  /  e)}[k)  denote  the  DFT  of  a(x/e)  on  the  grid;  that  is, 
Ma(x/e)](k)  =  f>(^)  e~2”lk'N • 

3=0  '  ' 


Then 


I- 


27V 


X 


t/2)\ 


27V— 1 


3=0 

N-l 


27 TJ 

2N -e/2 

2? rj 
Ne 


—2iri  —  i 


(2k)  =  Y. 

3=0  V“"' 

a(—^J  [e~2nijk/N  +  e-27T7(j+JV)fc/JVj 


Q  |  [e~27r ijk/N  e~2mjk/N e~2mk~^ 


3=0 

=  2  TN[a(x/e)}(k), 


so  that  the  even  coefficients  of  J-~2n 


e/2 


are  just  those  of  T n[o-{x  /  e)\.  Also, 


T- 


2N 


X 

JJ2 


27V- 1 


3=0  V  '  ' 


TV— 1 

E‘ 

3=0 

N-l 

l=o 

TV— 1 


2ttJ 

Ne 


z~2™^3  +  e-2wi^±l(j'+TV)' 


Ate 


1  +  e"27ri^ 


3=0 
=  0 


so  that  all  odd  coefficients  vanish.  These  equalities  give  2.15. 


□ 
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